@mapleus It appears that you just want to save results computed in each iteration in a loop (what you call a cycle).

Example:

restart;
for i from 1 to 5 do
  u:=(x1,x2,z)->x1+x2*z+9*z^i;
  one_proc:=proc(x1,x2) u(x1,x2,1) end proc;
  two_proc:=proc(x1,x2) u(x1,x2,2) end proc;
  res:=fsolve([one_proc,two_proc],[0,0]);
  print(res);
  C[i]:=subs([x1,x2]=~res,6.123*x1+4*x2)
end do:

eval(C);

@kle8309 Am I misunderstanding something:

You can answer your own question and then select it as the best answer?

Anyway this seems to be what happened here. Silly!

@mapleus You could do:

restart;
c:=1234567;
save c,"G:/MapleDiverse/test.m";
restart;
c; #c is just c
restart;
read "G:/MapleDiverse/test.m";
c;

@Markiyan Hirnyk You are right. This feature doesn't exist in Maple 15.

Where does the Plot2D come from? It should be plainly plot.

@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