Do you have any reason whatsoever to think that finding an analytical solution is possible?

Using phi and phi0 instead of your 2D versions:

eq1:=diff(phi(tau), tau)-1-a(tau)^2*cos(2*phi(tau))-a(tau)^2+cos(phi0-phi(tau))/a(tau)-cos(2*phi(tau)) = 0;
eq2:=diff(a(tau), tau)-a(tau)^3*sin(2*phi(tau))+a(tau)-sin(phi0-phi(tau))-a(tau)*sin(2*phi(tau)) = 0;
#dsolve({eq1,eq2}); #Doesn't seem to end
#Just try setting a=1 to see the complexity of this:
eq3:=eval(eq1,a=1);
dsolve(eq3);
value(%);

@mapleus The problem lies with fsolve, or rather with the difficulty of finding an accurate numerical solution to something like x^3 = 0. The derivative at the solution is zero. Thus Newton's method, which relies on the derivative at the approximate point, will get into trouble. Surely, fsolve is set up to handle that kind of thing, but apparently fails to find a solution that it itself is satisfied with.
Here is an example obtained by taking n:=3 and i:=1.

restart;
n:=3: i:=1:
u:=(x1,x2,z)->(x1+x2*z+9*z^i)^n;
p:=proc(x1,x2) option remember; local res;
  res:=[u(x1,x2,1),u(x1,x2,2)];
  #print([_passed,res]); It will print an awful lot if you remove #
  res
end proc:
p1:=proc(x1,x2) p(_passed)[1] end proc;
p2:=proc(x1,x2) p(_passed)[2] end proc;
fsolve([p1,p2],[0,0]); #Responds quickly unevaluated
T:=op(4,eval(p)): #The remember table of p (huge)
#eval(T); #To see it remove #
L:=entries(T,nolist):
AL:=map(x->abs(x[1])+abs(x[2]),[L]): #The sum of absolute values of the error in making zero
min(AL); #The minimal error
##Or you could sort AL
SAL:=sort(AL):
SAL[1..5];
nops(AL);
You see that fsolve is rather reluctant to declare success!
It is not too difficult to find the input corresponding to the minimal error (i.e. the best point (x1,x2) ), but since this is only an illustration I stop here.

@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.