@nm How about using showstat:

showstat(dsolve);

The upper limit found above can be computed using notation from the same wikipedia link.
I'm now using the original matrices.

restart;
A := Matrix([[-1, 0, 0, (1/2)*sqrt(2), 1, 0, 0, 0], [0, -1, 0, (1/2)*sqrt(2), 0, 0, 0, 0], [0, 0, -1, 0, 0, 0, 1/2, 0], [0, 0, 0, -(1/2)*sqrt(2), 0, -1, -1/2, 0], [0, 0, 0, 0, -1, 0, 0, 1], [0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, -(1/2)*sqrt(2), 0, 0, (1/2)*sqrt(3), 0], [0, 0, 0, 0, 0, 0, -(1/2)*sqrt(3), -1]]);

b := Vector([0, 0, 0, 0, 0, 10000, 0, 0]);
#Splitting A into A = L+D1+U:
D1:=Matrix(8,(i,j)->`if`(i=j,A[i,j],0));
L:=Matrix(8,(i,j)->`if`(i>j,A[i,j],0));
U:=Matrix(8,(i,j)->`if`(i<j,A[i,j],0));
##Convergence depends on the size of the absolute values of the eigenvalues of the following matrix M. They should all be less than 1.

M:=(D1+omega*L)^(-1).(omega*U+(omega-1)*D1);
ev:=LinearAlgebra:-Eigenvalues(M);
r:=solve(ev[1]-1,omega);
evalf(r); #Returns 1.136469738
# The issue is not settled yet, but from the 6 eigenvalues thar are equal to omega-1 it follows that 0<omega < 2 must be required. Further omega needs to be less than r since ev[1] increases with omega for omega>1. Thus 0<omega<r is required. But is it sufficient. We simply plot:

plots:-complexplot([ev[1],ev[2]],omega=0..r,scaling=constrained);
#The same animated:
plots:-animate(plots:-complexplot,[[ev[1],ev[2]],omega=0..s,scaling=constrained,thickness=3],s=0..r);



#We notice that the two first eigenvalues stay within the unit disk in the complex plane.

@Carl Love The subs trick is nice, but should be used with care:

restart:
A:=<1,5,1|1,3,1|6,2,1>;
subs(1=0,A); #OK

S:=series(exp(x),x);
subs(1=0,S); #Probably not expected

@mapleus What you need instead of unapply is to use output=listprocedue: You get a list of procedures (well actually equations whose right hand sides are procedures). In our case one of them will be u(z)=proc(...) ... end proc.
You want to get that procedure. Let us call it U (capital U to avoid confusion with the name u).

restart;
B:=1:
q:=1000:
l:=1:
n:=4.7:
V:=1:
M_F:=q*(2*l*(z-l)-z^2/2):
M_1:=piecewise((z<l), l-z, 0):
M_2:=2*l-z:
 

 one_proc:=proc(X_1,X_2) local res,M_finish,M_use,eq,cond,sol1,U;
 M_finish:=M_1*X_1+M_2*X_2+M_F:
 M_use:=signum(M_finish)*abs((M_finish)^n):
 eq := diff(u(z), `$`(z, 3)) = B/V*M_use:
 cond := u(0) = 0, (D(u))(0) = 0, ((D@@2)(u))(0) = 0;
 sol1 := dsolve({cond, eq}, numeric,output=listprocedure); ##
 U:=subs(sol1,u(z)); #Capital U
 res:=U(2*l); #Ditto
 print([_passed],res);
 res
 end proc;

 two_proc:=proc(X_1,X_2) local res,M_finish,M_use,eq,cond,sol1,U;
 M_finish:=M_1*X_1+M_2*X_2+M_F:
 M_use:=signum(M_finish)*abs((M_finish)^n):
  eq := diff(u(z), `$`(z, 3)) = B/V*M_use:
  cond := u(0) = 0, (D(u))(0) = 0, ((D@@2)(u))(0) = 0;
  sol1 := dsolve({cond, eq}, numeric,output=listprocedure); ##
  U:=subs(sol1,u(z)); #Capital U
  res:=U(2*l); #Ditto
  print([_passed],res);
  res
  end proc;

fsolve([one_proc,two_proc],[500,500]);
#OUTPUT:
[500, 637.17680634658297]

####### Time to wonder why X_1=500 ! Well, M_1 is zero for z = 2*l, so any value for X_1 is a solution together with the correct X_2.
#############

@Bryon This is part of what I see at

http://www.mapleprimes.com/questions/202870-Solving-Systems-Of-Nonlinear-Equations

@9colai I don't know what is the problem in your case, I would have to see the main document or at least the part involving the differential equation and any associated assignments to parameters appearing in that ode.

However, I can generate a similar error like this:

ode:=diff(y(t),t,t)+y(t)=0;
ICS:=Theta(0) = 2/3*Pi, D(Theta)(0) = 0;
dsolve({ode,ICS});

@9colai With a large right hand side my inclination would be to use numerical solution. However, if you for some reason want a formula for the solution and not just a numerical procedure, you can do as follows (where I compare with the numerical solution):
Assuming that your ode is simply named 'ode', and your initial conditions are ICS.
#For comparison:
resnum:=dsolve({ode,ICS},numeric);
plots:-odeplot(resnum,[t,Theta(t)],0..10); p1:=%:
#Now find a formula for general rhs equal to f(t)
res:=dsolve({lhs(ode)=f(t),ICS});
odetest(res,[lhs(ode)=f(t),ICS]); #OK
resHOM:=value(eval(rhs(res),f=0)); #The solution to the homogeneous equation
#Particular solutions with rhs = A*sin(a*t) and rhs= A*a*cos(a*t), respectively, and satisfying homogeneous initial conditions:
resS:=value(eval(rhs(res),f=(t->A*sin(a*t))))-resHOM: #rhs = A*sin(a*t)
resC:=value(eval(rhs(res),f=(t->A*cos(a*t))))-resHOM: #rhs= A*a*cos(a*t)
F:=rhs(ode):
indets(F,specfunc(anything,sin));
indets(F,specfunc(anything,cos));
#The next 3 lines are just an attempt at an explanation of what is going to follow:
match(op(3,F)=A*sin(a*t),t,'s');
s;
evalf(eval(resS,s));
#Now we find the solution as a sum of the solution to the homogeneous ode with ICS and a sequence of particular solutions corresponding to different right hand sides, but each satisfying homogeneous initial conditions. We go through all the terms in F:
sol:=evalf(resHOM):
for i in F do
  if match(i=A*sin(a*t),t,'s') then
    sol:=sol+evalf(eval(resS,s))
  elif match(i=A*cos(a*t),t,'s') then
    sol:=sol+evalf(eval(resC,s))
  else
    error "no match"  #Shouldn't happen
  end if
end do;
plot(sol,t=0..10); p2:=%:
plots:-display(p1,p2);

@9colai Copying and pasting your expression and assigning it to to u I see a factor 10^698, which is huge.
u has 17 terms. All but one are innocuous. The first has the ominous  10^698.
To display the terms:
for i to 17 do op(i,u) end do;
#Plotting
plot(add(op(i,u),i=2..17),t=0..60); #No problem
plot(op(1,u),t=0..60); #Problem

Solution:

w:=expand(u);
plot(w,t=0..60):

Except that now a=1 this is the same question you asked in

http://www.mapleprimes.com/questions/202911-Warning-Solutions-May-Have-Been-Lost

@mapleus Notice that although certainly a^n1*a^n2=a^(n1+n2) the identity you like to use is
a^(n1+n2) = a^n1*abs(a)^n2, which is not correct, but amounts to a reinterpretation of which integral you want to make zero. Also if n1 is even the problem is as before. Thus you should pick n1 odd to make sure the sign can change in the interval. But then a cleaner (but equivalent) interpretation for F(z) is
F(z)^n = signum(F(z))*abs(F(z))^n.

@Kitonum Thanks, you caught a special case: If h2=h1+1 you get a division by zero.
It seems that that can be handled by limits:

TresA:=limit(rhs(Tres2),h2=h1+1);
SresA:=limit(rhs(Sres2),h2=h1+1);
odetest({T(y)=TresA,sigma(y)=SresA},{Eq2,Eq3} union eval({bcs2},h2=h1+1));

@9009134 You can obtain the six solutions like this:

1. Pick as required {f2(0) = 0, f3(0) = 0, (D(f1))(0) = 0}. (Required by the differentiation resulting in ode).
Add to that one (only one) of {f2(1) = 0, f3(1) = 0, (D(f1))(1) = 0}. It doen't matter which since the remaining two will be satisfied because of conditions R1 and R2 (se above).
You now need to add 3 conditions from the available rest, i.e. those not having any of the conditions so far mentioned. Since there are 4 remaining conditions, you get 4 solutions that way.
2. Now still picking the required {f2(0) = 0, f3(0) = 0, (D(f1))(0) = 0} you add the 4 remaining conditions mentioned in point 1, i.e. this time you don't take any of {f2(1) = 0, f3(1) = 0, (D(f1))(1) = 0}. That gives you one solution. Thus so far you have 5 solutions.
3. Finally pick as required {f2(1) = 0, f3(1) = 0, (D(f1))(1) = 0}. If you pick any of {f2(0) = 0, f3(0) = 0, (D(f1))(0) = 0} you will have them all because of R1 and R2, so you won't be getting any new solutions. They will be the same as the 4 solutions from point 1.
So to get the last solution add to the required {f2(1) = 0, f3(1) = 0, (D(f1))(1) = 0} the 4 conditions left after removing {f2(0) = 0, f3(0) = 0, (D(f1))(0) = 0}.

You got 6 solutions.

@9009134 

Each of the following 13 sets

[{f1(0) = 0, f1(1) = 0, f2(1) = 0, (D(f3))(1) = 0}, {f1(1) = 0, f2(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f2(1) = 0, (D(f3))(0) = 0}, {f1(0) = 0, (D(f1))(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(1) = 0, (D(f1))(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f3(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f3(1) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, (D(f1))(1) = 0, (D(f3))(0) = 0}, {f1(0) = 0, f1(1) = 0, (D(f1))(1) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f3(1) = 0, (D(f3))(0) = 0}, {f1(1) = 0, f3(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f2(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}]

union

{f2(0) = 0, f3(0) = 0, (D(f1))(0) = 0}

and each of the following

[{f1(0) = 0, f2(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(1) = 0, f2(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(1) = 0, f3(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, (D(f1))(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f3(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(1) = 0, (D(f1))(0) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f2(0) = 0, (D(f3))(0) = 0}, {f1(0) = 0, f1(1) = 0, (D(f1))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, (D(f1))(0) = 0, (D(f3))(0) = 0}, {f1(0) = 0, f1(1) = 0, f2(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f3(0) = 0, (D(f3))(1) = 0}, {f1(0) = 0, f1(1) = 0, f3(0) = 0, (D(f3))(0) = 0}]

union

{f2(1) = 0, f3(1) = 0, (D(f1))(1) = 0}

There actually only 6 different solutions. The reason is that these two relations:
R1:=-(D(f1))(0)+(D(f1))(1) = -(1/7)*f3(0)+(1/7)*f3(1)
R2:=-f2(0)+f2(1) = (11/7)*f3(0)-(11/7)*f3(1)
can found from
solve({sys[1],ode},{diff(f1(x),x,x),diff(f2(x),x)});
and then:
R1:=map(int,%[1],x=0..1,continuous);
R2:=map(int,%%[2],x=0..1,continuous);

Could you give us an explicit example?

To find solutions for which x>0,y>0 I believe it will be necessary to use numerical solution, but then of course you need all the parameters a, q1, ... to be concrete.