@Preben Alsholm You said, my problem can by a change of variables y(x) = A*v(z), x=L*z, then how to transfer the second order differential y''(x) to v''(z)?  because I don't know how do you calculate B= 0.549269696597376e-5 in the toy problem?  

you have mentioned the ode:

ode:=diff(y(x),x,x)=-abs(y(x))^(1/3)+x;

Can you solve it by code in maple?  Is it using the simple shooting method?  I have a lot of places confused.  Thank you.

@Preben Alsholm Hey, Preben.  I appreciated that you left too many words for me, but I can't find the answer what I need.  Sorry about that...  The simple shooting method is mentioned in other paper, it can be approved to solve this differential equation.  But I can't use it in maple.  I really don't  find the answer for it :

f''(x)=piecewise(f(x)>=0, 3.019*(10^(-5))*f(x)^0.337, [-3.019*(10^(-5))*(-f(x))^0.337])+9.541*(10^(-13))*x

f(0)=0;  f'(2945)=0        0<x<2945 is very important.

Thank you for your any help in advance.  Anyway, If you are also willing to help me...

@Preben Alsholm  I am waiting for your perfect solution.  Absolutely, you can figure it out.  I trust you can give the help what I need.  And I will spend more time on understanding the code for the toy problem,  Wish I may find something new.  In one word, thank you so much!   Good luck to us.

@Preben Alsholm  Hey, Preben.  Probably I didn't  clarify this problem. So I explain the equation what I focus on in order to find the solution.   

Firstly, I derived the governing equation shown below:

y''(x)-3.019*10^(-5)*u=0, 

where,   u=y(x)^0.337 come from the experience data.  There isn't part for the negative coordinate. 

After adding the new term [9.541*10^(-13)]*x on the right-hand side of ode such as y''(x)-3.019*10^(-5)*u=9.541*10^(-13)*x, we have to use the negative part of y(x). So I change the graph of u like this:

Due to the piecewise expression isn't differentiable at zero, so I revised this ode again. At the same time, we definitely make sure y(x) negative. I made u like this, it is available and suitable for my research. 

so, the final revised ode is shown below,

y''(x)-3.019*10^(-5)*u=9.541*10^(-13)*x, 

where,  u=-[(-y(x))^0.337] .   The boundary condition, y(0)=0, y'(2945)=0.

Can you help me solve this ode?    Thank you for your help in advance!!!

@Preben Alsholm In fact, the revised code I run like this in the below:

restart;

Digits := 15;

u := piecewise(y(x) < 0, -(-y(x))^0.337, 0);

ode2 := diff(y(x), x, x)-0.3019e-4*u = 9.542*10^(-13)*x;

bcs := y(0) = 0, (D(y))(2945) = 0; res2 := dsolve({bcs, ode2}, numeric, method = bvp[midrich], approxsoln = [y(x) = (1/2945)*x^2-2*x]);

plots:-odeplot(res2, [x, y(x)])

Using the u=piecewise(y(x) < 0, -(-y(x))^0.337, 0) , isn't u=piecewise(y(x)<0,0,-(-y(x))^0.337).  That was a slip of the pen.  I apologize for my mistake made you confusion. 

Can you help me solve it by maple?  I really appreciate your helpful reply every time! 

@Preben Alsholm  Yes. y(x) can change sign, this piecewise is continuous, but not differentiable at zero.

piecewise(y(x)>=0,y(x)^0.337,-((-y(x))^0.337))

So can we revise the ode in order to solve the equation?  Now I knew y(x)<0,  so the new ode is shown below:

y''(x)+0.00003019*[-y(x))]^0.337=[9.541*10^(-13)]*x 

The boundary condition always keeps same as before.  

Because we knew y(x) absolutely was negative, so I can change the sign of the term [-y(x)]^0.337 instead of y(x)^0.337.   

I tried to solve it in maple using the code written by you but showed me an error. 

restart;  

Digits:=15:  

u:=piecewise(y(x)<0,0,-(-y(x))^(0.337)));  

ode2:=diff(y(x),x,x)-0.00003019*u=9.542*10^(-13)*x;  

bcs:=y(0)=0,D(y)(2945)=0;  

res2:=dsolve({ode2,bcs},numeric,method=bvp[midrich],approxsoln=[y(x)= (1/2945)*x^2-2*x]);  

plots:-odeplot(res2,[x,y(x)]);

Would you like to give me some advice about it? Thank you for helping me all the time and your patience.  Sincerely! 

@Preben Alsholm  Hey, Preben.  I thought about the problem for couple days, I appreciate that you give me the evidence which proves y(x)<0. In the beginning, I considered the graph matched what I predicted. However, the equation is different from the initial equation (y''(x)-a*y(x)^0.337=b*x).  Because the graph is plotted by y''(x)=b*x.  At the same time, I am pretty sure that the boundary condition is correct. (y(0)=0, y'(2945)=0).  The governing equation for the model of the reinforced concrete simple beam which addressed the nature of the bond between a tendon and concrete in the transfer zone. Where y(x) is the slip of the beam, x=0 at the center of the beam, x=2945 at the end of the beam.   According to the situation of y(x),  I changed the governing equation such as:

y''(x)-3.019*10^(-5)*piecewise(y(x)>=0,y(x)^0.337,-((-y(x))^0.337))=9.541*10^(-13)*x

Boundary condition:

y(0)=0, y'(2945)=0 

Can you help me solve this changed governing equation in Maple? Thank you for your any help in advance!!!

@Kitonum GOOD JOB!  Thank you for your answer very much, because it gives me more thinkings and ideas. Thanks.

@Preben Alsholm 

Hey, Preben. First of all, I appreciate that you give me you 100% support for every time. You really don't how much important of your reply to me!  I have one question about your code.  Why do you set "approxsoln=[y(x)=tanh(x)]", I used to find the answer in Maple Online Help which means assuming a starting point for the problem.  But I still don't know it totally.  Can I make y(x)=0, or anything else? Thank you for your answer one more time!  Sincerely.

@JohnS The physical problem is the bond-slip relationship at the tendon and concrete interface, the detail is shown in the comment below. 

@Rouben Rostamian  The picture shown below is the coordinate system of the physical model of the reinforced concrete. simple beam.

@Rouben Rostamian  Thank you for your perfect answer.  This is a model of the reinforced concrete simple beam which addresses the nature of the bond between a Tendon and concrete in the transfer zone.  Where y(x) is the slip at point x, slip y(x) always happen at the tendon and concrete interface.  Setting the initial point (x=0) at the center of the beam, x=L at the end of the beam. absolutely, y(0)=0 and y'(0)=0.

I appreciate that your answer can help my work go further. Thank you so much! 

@Preben Alsholm  Thank you for helping me! You are pretty nice to me,  for giving one more answers to me.  I appreciate that you focus on this kind of question continually.   You offer me the new method of numerical solution in maple. That is awesome!

@Rouben Rostamian  Your answer is the best answer that I have ever seen. Thank you for helping me a lot!!!  

Can I know where do you find the elementary function for solving this differential equation, is it very famous?  Because I never see it before.  By the way, my major is Civil Engineering.  

And I saw your note when L < L__crit, Maple's numerical ODE solver can figure it out. Can you show me the code of numerical solution?  Actually, The closed form is what I want, But the results y(x<L-L__crit)=0 in the exact solution confused me, I don't think it is equal to 0.  So I need the numerical solutions for 0<x<L-L__crit if there is no exact solution for this interval.  Hopefully, you can help me find the solution for the remaining case(0<=x<=L-L__crit).  Thank you in advance.

Thank you for your time and your assistance in this matter! 

Sincerely,

oceanxinlei