@Preben Alsholm 
Yes, I again made a mistake. Thank you!
it was necessary to take into account that sin(x)<=1.
then the code to change: 
instead
r_min:=0.3-step:
must be
r_min:=0.152-step:
still wondering how to get the value of the integrated functionat a specific point.

@Preben Alsholm 
you are right
thanks for the correction
one iterationworks too (notgiving errors) ifI use an olderversion of their program(without your patchaboutsubs - r_ravn:=subs(F,r(s));).

Tell me, please,how do I"get"the solution?
Isolve the system, resulting in the solution forr(s),buthow to get the valueof this solutionat a specific point,for example, r(s_end)???

the initial error iscorrected:

restart:
R:=0.3:
theta_min:=Pi/6:
theta_max:=Pi/2:
betta_max:=evalf(Pi/180*80);
p:=2*10^5:

theta0:=s->Pi/3/s_end*s+Pi/6:
r0:=s->R*sin(theta0(s)):
s_end:=evalf(R*(theta_max-theta_min)):

sol1:=solve({sin(betta_max)=c/r0(0)},{c});
const1:=0.1477211630;

betta0:=s->arcsin(const1/r0(s)):
betta:=s->arcsin(r(s)/r0(s)*sin(betta0(s))):
A:=s->cos(betta(s))/cos(betta0(s)):
T1:=s->rT1(s)/r(s):
T2:=s->T1(s)*tan(betta(s))^2:

step:=0.01:
delta:=0.001:
for i from 1 to 1 do
r_min:=0.3-step:
rT1_n:=p*Pi*r_min^2/2/Pi/sin(theta_min):

sys := diff(rT1(s), s)-A(s)*T2(s)*cos(theta(s)),diff(theta(s), s)-A(s)/T1(s)*(p-T2(s)*sin(theta(s)/r(s))),diff(r(s),s)-A(s)*cos(theta(s)),diff(z(s),s)-A(s)*sin(theta(s));
fcns := {rT1(s),theta(s),r(s), z(s)};
F := dsolve({sys,rT1(0)=rT1_n, theta(0)=theta_min,r(0) = r_min, z(0) = 0}, fcns, numeric,output=listprocedure):
r_ravn:=s->subs(F,r(s)):
if abs(evalf(r_ravn(s_end))-R)=delta then break:
print(r_min):
end if:
end do:

r_ravn:=s->subs(F,r(s));

evalf(r_ravn(s_end));
evalf(r_ravn(0));
plot([r_ravn(s),r(s)],s=0..s_end);

Carl Любовь 
Thanks for the help!
signum I did not finish by mistake: (
*signum(M_x(z))

But I can't see the result and graph, for example, I write:
evalf(X), the answer is : X
I write :
sol_4:=fsolve(sol_3);
M_xx:=z->subs(sol_4,M_x(z)):
plot(M_xx(z),z=0..3*l);
answer:
Error, (in M_xx) invalid input: subs received fsolve(sol_3), which is not valid for its 1st argument
Or I should go for the procedure?

@Carl Love 

thanks!

 restart;
 with(DirectSearch):
 B:=1: 
 q:=1: 
 l:=1: 
 n:=4.7:
 V:=0.05:
 N:=1200:
 
 
 kappa:=Vector(N+1,[]):
 theta:=Vector(N+1,[]):
 u:=Vector(N,[]):
 M:=Vector(N,[]):
 Z:=Vector(N,[]):
 
 M_F:=q*(6*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_m:=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_finish:=unapply(M_m,z):
 
 kappa_old:=0:
 theta_old:=0:
 u_old:=0:
 M_old:=0:
 
 
 step:=6*l/N:
 u[1]:=0:
 kappa[1]:=0:
 theta[1]:=0:
 
 
 
 
 for i from 2 to N do
 
 z:=i*step:
 kappa_new:=kappa_old+B/V*(M_finish(z))^n*step:
 
 theta_new:=theta_old+1/2*(kappa_old+kappa_new)*step:
 
 u_new:=u_old+1/2*(theta_old+theta_new)*step:
 
 Z[i]:=z:
 kappa[i]:=kappa_new:
 theta[i]:=theta_new:
 u[i]:=u_new:
 kappa_old:=kappa_new:
 theta_old:=theta_new:
 u_old:=u_new:
 
 end do:
 
 
 u[N/6]=0;
 u[N/3]=0;
 u[N/2]=0;
 u[2*N/3]=0;
 u[5*N/6]=0;
 u[N]=0;
 

Please advise how to solve the resulting system.

@Preben Alsholm 
Thank you! Yes, in a loop.

If possible, I would like to ask another question, about the solution of a system of nonlinear equations.
If in your example, the equation to transform in nonlinear

u:=(x1,x2,z)->(x1+x2*z+9*z^i)^(n); where n is float

That solution does not work, even if "n" is an integer.
In your first remark (about the integrals), the system had a solution. This is similar, but no solution. Perhaps I don't understand the internal mechanism.
But I really want to learn how to solve systems of nonlinear equations with float indicator, such as 

u:=(x1,x2,z)->(x1+x2*z+9*z^2)^(n)

@Preben Alsholm 

Thank you!
but using the cycle
I won't be able to do c:=1234567

Does this mean that the procedure that defines the values X_1 and X_2 can not be used in a cycle?
It is possible to arrange something like that (any example):
(Yes, this does not apply to the original question, but I'd like to understand whether it is possible and how to use procedures in cycles.)

k:=1:
n:=4.7: 

for i from 1 to 100 do

u:=k*(F(X_1,X_2))^n:
#where F is any function of X_1 and X_2, and u contains a "variable" z = l..2 * l.

one_proc:=proc(X_1,X_2) local res;
  
  res:=evalf(u(l));
  print([_passed],res);
  res 
end proc; 

two_proc:=proc(X_1,X_2) local res;
  
  res:=evalf(u(2*l));
  print([_passed],res);
  res 
end proc; 

#???

c:=a*X_1+b*X_2:

C[i]:=c:
k:=k+1:

end do:

@Preben Alsholm 

I have:

eq := diff(u(z), `$`(z, 3)) = B/V*M_use:
cond := u(0) = 0, (D(u))(0) = 0, ((D@@2)(u))(0) = 0;

I have differential equationof the third order.
The function u I can be represented asa triple integral.And I wantto consistently calculateone from the other, smashing the siteof integration on intervals and use cycle.

But I don't know how in the cycle to go to six procedure sand how to use the result obtained with their help.

@Preben Alsholm 

In the beginning you said:

restart; 
B:=1: 
q:=1*10^3: 
l:=1: 
n:=4.7: 
M_F:=2*l*(z-l)-z^2/2; #No need to make the Ms functions. Divided by q
M_1:=piecewise(z<l, 2*l*(z-l)+(l-z)/q-z^2/2, 2*l*(z-l)-z^2/2); #Divided by q
M_2:=2*l*(z-l)+(2*l-z)/q-z^2/2; #Divided by q
#one_ and two_int should be functions of X_1 and X_2, not z.
one_int:=proc(X_1,X_2) local res;
  res:=evalf(Int(B*(M_F+X_1*M_1+X_2*M_2)^n*M_1,z=0..2*l));
  print([_passed],res);
  res 
end proc; 
two_int:=proc(X_1,X_2) local res;
  res:=evalf(Int(B*(M_F+X_1*M_1+X_2*M_2)^n*M_2,z=0..2*l));
  print([_passed],res);
  res 
end proc;
#Printing in these procedures lets us know what is going on.
##The following prints out for quite a while and then stops printing, but it keeps evaluating. I interrupted the computation. 
fsolve([one_int,two_int],[-20,19]);

Please answer the question:
Is it possible to make so that the results are automatically usedi n further calculations?

For example, I have :   c:=a*X_1+b*X_2
andafter restarting thecodeI havewasthe result с="the number"

Preben Alsholm Yes, the original question was about solving the system of nonlinear equations. Now I want to try to solve the same problem in a different way. However, the problem of solving the system of nonlinear equations remains. Only in cycle. With the method of use of the procedures I met thanks to you. But my experience with maple is very small and in some cases can not do without help.

Preben Alsholm 
Can you help, please?
Then anothererror "ErrorVectorindex out ofrange"
Unfortunately, I don't knowhow to fix it.

@Preben Alsholm

 Тhank you very much!


Very well-written!

But, unfortunately, the conditions laid down in the solution U(k*l))=0 or U(k*l))=0.001 (k=1..6) are satisfied only for k=2 (the system of 2 equations) and almost true for k=3 (the system of three equations).

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.

In cycle
at each step, I'll have u(X_1,X_2,...,X_6). And conditions u(l)=0, u(2*l)=0,..,u(6*l)=0.
But I have a problem.
I don't know how i can  apply in the cycle  the procedure search X_1,...X_6 at each step (n=4.7) and transfer these values to further account on the same step.

restart;
 B:=1:
 q:=1000:
 l:=1:
 n:=4.7:
 V:=0.05:

 kappa:=Vector(N,[]):
 theta:=Vector(N,[]):
 u:=Vector(N,[]):
 M:=Vector(N,[]):
 Z:=Vector(N,[]):

 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:

 kappa_old:=0:
 theta_old:=0:
 u_old:=0:
 M_old:=0:


N:=1000:
step:=6*l/N:

 for i from 1 to N do
 kappa_new:=kappa_old+B/V*(M_finish)^n*step:
 theta_new:=thena_old+1/2*(kappa_old+kappa_new)*step:
 u_new:=u_old+1/2*(theta_old+theta_new)*step:


##If you use 6 procedures, it is necessary to combine their solutions into one "sol".
#Or can all be combined into a single procedure?
#How would the meaning should be implemented like that. Within the meaning of. How to write real?
 sol:=fsolve({u_new(l)=0,u_new(2*l)=0,u_new(3*l)=0,u_new(4*l)=0,u_new(5*l)=0,u_new(6*l)=0},X_1,X_2,X_3,  X_4,X_5,X_6}):



 u_temp:=subs(sol,u_new):
 u[i]:=temp(i*step):
 M_temp:=subs(sol,M_finih):
 M_new:=1/2*(M_old+M_temp(i*step)):
 M[i]:=M_new:
 Z[i]:=i*step:
 kappa_old:=kappa_new:
 theta_old:=theta_new:
 u_old:=u_new:
 M_old:=M_new:

 end do:

u_plot := [seq([Z[i], u[i]], i = 1..N)]:
plot({u_plot},color=blue,thickness=2);

M_plot := [seq([Z[i], M[i]], i = 1..N)]:
plot({M_plot},color=green,thickness=2);

Preben Alsholm 
Тhank you very much!
in proxy time you really helped!
But I have  problems again.
My system for the general case contains 6 equations and 6 unknowns.
X_1, X_1,...,X_6
To find their values from the conditions:
res:=U(n*l) # =0 or 0.001
for examlpe 
res:=U(l),
res:=U(3*l) ...
But I have error.
Can suggest how to fix the error?

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:

one_proc:=proc(X_1,X_2,X_3,X_4,X_5,X_6) local res,M_finish,M_use,eq,cond,sol1,U;
 M_finish:=M_1*X_1+M_2*X_2+M_3*X_3+M_4*X_4+M_5*X5+M_6*X_6+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)); 
 res:=U(l); #Ditto
 print([_passed],res);
 res
 end proc;

two_proc:=proc(X_1,X_2,X_3,X_4,X_5,X_6) local res,M_finish,M_use,eq,cond,sol1,U;
 M_finish:=M_1*X_1+M_2*X_2+M_3*X_3+M_4*X_4+M_5*X5+M_6*X_6+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)); 
 res:=U(2*l); #Ditto
 print([_passed],res);
 res
 end proc;

three_proc:=proc(X_1,X_2,X_3,X_4,X_5,X_6) local res,M_finish,M_use,eq,cond,sol1,U;
 M_finish:=M_1*X_1+M_2*X_2+M_3*X_3+M_4*X_4+M_5*X5+M_6*X_6+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(3*l); #Ditto
 print([_passed],res);
 res-0.001
 end proc;

four_proc:=proc(X_1,X_2,X_3,X_4,X_5,X_6) local res,M_finish,M_use,eq,cond,sol1,U;
 M_finish:=M_1*X_1+M_2*X_2+M_3*X_3+M_4*X_4+M_5*X5+M_6*X_6+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(4*l); #Ditto
 print([_passed],res);
 res-0.001
 end proc;

five_proc:=proc(X_1,X_2,X_3,X_4,X_5,X_6) local res,M_finish,M_use,eq,cond,sol1,U;
 M_finish:=M_1*X_1+M_2*X_2+M_3*X_3+M_4*X_4+M_5*X5+M_6*X_6+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(5*l); #Ditto
 print([_passed],res);
 res-0.001
 end proc;

six_proc:=proc(X_1,X_2,X_3,X_4,X_5,X_6) local res,M_finish,M_use,eq,cond,sol1,U;
 M_finish:=M_1*X_1+M_2*X_2+M_3*X_3+M_4*X_4+M_5*X5+M_6*X_6+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(6*l); #Ditto
 print([_passed],res);
 res
 end proc;

fsolve([one_proc,two_proc,three_proc,four_proc,five_proc,six_proc],[100,100,100,100,100,100]);

@Preben Alsholm 
Thanks for the greatidea.
The solution was obtained, but there has been little changein the condition and my attempts decisions result in an error.
Please, help to correct errors.

 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_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);
 u:=unapply(sol1,z);
 res:=evalf(u(l));
 print([_passed],res);
 res
 end proc;



 two_proc:=proc(X_1,X_2) local res;
 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);
  u:=unapply(sol1,z);
  res:=evalf(u(2*l));
  print([_passed],res);
  res
  end proc;



fsolve([one_proc,two_proc],[500,500]):

Error, (in fsolve/sysnewton) procedures don't evaluate to numeric types.

@Preben Alsholm 
thanks for the detailedexplanation!!!



I have
A=Int(F(z)^n*h(z),z=0..2),
F(z):=(M_F+X_1*M_1+X_2*M_2), where M_F<0  if 0<z<2
And I will get the solution if  n is an integer
or n is float if F(z)>0 , z:=[0..2].

But I can combine these 2 conditions:
fot example, n - float.

F(z)^n = [F(z)^(n1)]*[abs(F(z))^(n2)]

{a^n1*a^n2=a^(n1+n2)}

where n1 is integer, n1 is float.

Can I do this?

Preben Alsholm 
Thanks!

I slightly changed initial conditions

B:=.1470494889e-12:

M_1:=piecewise((z<l), l-z,  0);

M_2:=2*l-z;

And now, when  n non-integer, for example, n=4.7 integrals take complex values. What does this mean?
Тhe actual solution does not exist,if to solve technical applied task?

For exampl:
[4.8019653063270831, -4.1602326280958588], -45.30522058416861298+62.35728651500813821*I

If n is an integer, then the result is a real number.