@mapleus When you switched from the original problem involving two integrals over z=0..2*l to the problem with
 eq := diff(u(z), `$`(z, 3)) = B/V*M_use:
 cond := u(0) = 0, (D(u))(0) = 0, ((D@@2)(u))(0) = 0;

which really has as solution a triple integral, it seems to me that you changed the problem quite a bit.

When you say:

I want to solve it "manually".
To divide the interval of integration on small segments (N=1000) and to calculate the integral on each of them, for example, by the method of rectangles or trapezoids.

then which integral are you talking about?

@mapleus Honestly I'm not at all sure what you are really trying to do. Are you actually still talking about the problem in the original question?

How come you don't post your equations using Maple syntax?

@digerdiga You may want to look at the answer by Robert Israel in the link

http://www.mapleprimes.com/questions/87735-Animation-Of-Nonlinear-Waves-With-Maple#answer87979

Notice that the  TWSolution found by PDEtools:-TWSolution cannot be made to satisfy your boundary condition u(0, t) = u(2, t).

@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