@digerdiga I tried

pds:-animate(t = .01,frames=100);

i.e. animating over a much shorter time interval.

What is the solution supposed to do?

Could you provide some details, e.g. an uploaded worksheet. I don't understand the problem.

The syntax x:=A[1,5}  would be wrong no matter what A is, but it is probably a typographical error(?).

@yangsong525525 Pexact is a plot of the exact solution. The symbol used is a cross. In the last line this plot is displayed together with p[2] and p[3].

To see the error you could continue with:

xval:=map2(op,1,data); #The x-values used
E:=evalf(eval~(rhs(exact),X=~xval)); #The corresponding Y-values of the exact solution
RK:=map2(op,2,data); #The Y-values of the RK solution
RK-E; #The list of errors
ERR:=`[]`~(xval,RK-E); #Same but x-values are included for plotting purposes
plot(ERR,style=point,caption="RK minus exact");

I tested your code in Maple 12.02: No problem. Similarly, no problem in Maple 18.02.

@mohitgupta1993 Although the print on the Maple screen of your sys may look good to you, pdsolve will not understand it. Use diff as I said.
Like this:

pde1:=-A*kappa3-(2*G-A)*diff(W(x,y),x,x)-2*G*(diff(W(x,y),y,y)+diff(Z(x,y),x,y))+A*diff(Z(x,y),x,y)=0;
pde2:=(A-4*G)*diff(Z(x,y),y,y)+(A-2*G)*diff(W(x,y),x,y)-2*G*diff(Z(x,y),x,x)=0;
sys:={pde1,pde2};
#WARNING: If you try the following commented line the computation may not stop.
#pdsolve(sys);
##The outcome of this special case should discourage most people:
eval(sys,{A=0,G=1});
pdsolve(%);

@mapleus The error was that X5 inside the procedures should be X_5.
But you are solving the very same differential equation 6 times for each sequence of X's.
This is costing a lot of time.
Here is a much faster way using the parameters option to dsolve.
The procedure p sets the parameters in the statement sol(parameters=[X_1,X_2,X_3,X_4,X_5,X_6]);
It outputs the list [U(l),U(2*l),U(3*l)-0.001,U(4*l)-0.001,U(5*l)-0.001,U(6*l)].
The procedure p remembers which input got which result (option remember).
That has been done so that the procedures Q[i], i=1..6, which all call p, don't make p compute more than just once (and not 6 times).
That was the good news.
The bad news is that fsolve didn't finish with a result, but returned unevaluated.

If you want printout then set infolevel[p]:=2 (notice 'userinfo' in p).

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:=piecewise((z<2*l), 2*l-z, 0):
M_3:=piecewise((z<3*l), 3*l-z, 0):
M_4:=piecewise((z<4*l), 4*l-z, 0):
M_5:=piecewise((z<5*l), 5*l-z, 0):
M_6:=6*l-z:
M_finish:=M_1*X_1+M_2*X_2+M_3*X_3+M_4*X_4+M_5*X_5+M_6*X_6+M_F:
M_use:=signum(M_finish)*abs((M_finish)^n):
eq := diff(u(z), z,z,z) = B/V*M_use;
cond := u(0) = 0, D(u)(0) = 0, (D@@2)(u)(0) = 0;
sol := dsolve({cond, eq}, numeric,output=listprocedure,parameters=[X_1,X_2,X_3,X_4,X_5,X_6]);
U:=subs(sol,u(z));

p:=proc(X_1,X_2,X_3,X_4,X_5,X_6) option remember; local res;
   sol(parameters=[X_1,X_2,X_3,X_4,X_5,X_6]);
   res:=[U(l),U(2*l),U(3*l)-0.001,U(4*l)-0.001,U(5*l)-0.001,U(6*l)];
   userinfo(2,p,[_passed],evalf[3](res));
   res
 end proc:
for i from 1 to 6 do
  Q[i]:=subs(ii=i,(proc(X_1,X_2,X_3,X_4,X_5,X_6) p(_passed)[ii] end proc))
end do;
p(100,100,100,100,100,100);
infolevel[p]:=2:
Q[3](110,100,100,100,100,100);
infolevel[p]:=1: #Now no intermediate results, but you could try keeping it at 2.
fsolve([seq(Q[i],i=1..6)],[100,100,100,100,100,100]);
#I also tried with no great result:
Optimization:-LSSolve([seq(Q[i],i=1..6)],initialpoint=[100,100,100,100,100,100]);
##
Certainly the starting point [100,100,100,100,100,100] is not anywhere near a solution:
sol(parameters=[100$6]); #Setting parameters to [100,100,100,100,100,100]
plots:-odeplot(sol,[z,ln(abs(u(z)))],0..6*l); #Notice the logarithm

@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