@nm I think you are being rather unreasonable.
When I do
?pi
it brings up a help page, which starts out saying
"For information on 3.14159..., see Pi."

I think that that is as helpful as it gets!

And that
?foo
brings up a help page about Units,length, (foot)
just shows that the help system is trying to be helpful.

In your worksheet you don't attempt any solution of the system (sys), but you are finding equilibria by using solve. I don't get any error at all, but the output:



i.e. infinitely many equilibria, here parametrized by v, x, and z.

If you want to solve sys, you will undoubtably have to do it numerically. For that the constants m and q have to be made concrete and initial (or boundary) values must be known.

@sharena2 Why do you think results will get better using shooting?

However, if you insist, I suggest shooting from eta=L, which will only require two parameters.
For comparison I continue from the code I gave above. Thus 'sol' is the output from dsolve/bvp.

#Shooting from the right with the 2 parameters fL and vL:
ICSL:={F(L) = 0, H(L) = 2, f(L) = fL, g(L) =-fL, u(L) = 0,v(L)=vL};
resL:=dsolve(ODE union ICSL,numeric,parameters=[fL,vL],output=listprocedure);
#We shall need these two:
ff,uu:=op(subs(resL,[f,u](eta)));
#Getting parameter values from sol (no shooting):
pp:=subs(sol(L),[f(eta),v(eta)]);
#Quick graphical test:
resL(parameters=pp):
plots:-display(Array([[seq(plots:-odeplot(resL,A[i],0..L),i=1..3)],[seq(plots:-odeplot(resL,A[i],0..L),i=4..6)]]));

#There are 2 parameters  (fL and vL) and 2 requirements: f(0) = 0 and u(0) = 1
pL:=proc(fL,vL) option remember;
  resL(parameters=[fL,vL]);
  [ff(0),uu(0)-1]
end proc:
p1:=proc(fL,vL) pL(_passed)[1] end proc:
p2:=proc(fL,vL) pL(_passed)[2] end proc:
#Testing pL on the parameters from sol (pp):
pL(op(pp)); #Should be close to zero
ppL:=fsolve([p1,p2],pp);
pp;
ppL-pp;
#With no so good a guess:
fsolve([p1,p2],[.4,0]);
resL(parameters=%):
plots:-display(Array([[seq(plots:-odeplot(resL,A[i],0..L),i=1..3)],[seq(plots:-odeplot(sol,A[i],0..L),i=4..6)]]));

You normally would do something like (example)

if not type(aaa,numeric) then error "%1 should be of type numeric",aaa end if;

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