@oceanxinlie Your last  question:

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

Yes, a simple shooting method is used and I already gave the full code in Maple in the reply with the title "Not quite yet" at the end. Notice that the code starts with a restart command. That is a clear indication that everything is in there, nothing is left out except a detailed explanation of the code. I didn't care to do that since I hadn't (and still haven't) been successful in using it on your system transformed to a similar form. That transformation I shall give in a second reply.

@oceanxinlie The ode given by

ode:=diff(y(x),x,x)=abs(y(x))^0.337*signum(y(x))+b*x;

  is the same as the one you mention (after the transformation I mentioned above)

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

because
abs(y(x))^0.337*signum(y(x)) = y(x)^0.337 for y(x) >= 0
and
abs(y(x))^0.337*signum(y(x)) = -(-y(x))^0.337 for y(x) < 0.

@oceanxinlie Solving the problem

ode:=diff(y(x),x,x)=abs(y(x))^0.337*signum(y(x))+b*x;

considered on x = 0..1 using the shooting method shown above begins to run into problems already at values of b less than 0.1.
The b you have (once your ode is transformed as mentioned above) is 0.549269696597376e-5. For that value of b the graph of  y(0) as a function of y(1) shows what looks like a discontinuity where the sign changes:

It is probably not really a discontinuity, but "only" a numerical problem. It is not surprising, however, that the derivative of y(0) w.r.t. to the parameter y(0)=y0 seems to be infinite at a point where y = 0 because the exponent of y(x) in the ode is 0.337 < 1.
Surely it doesn't help that at x = 1 the value of y(1) = y0 must be chosen close to zero at about -2.5*10^(-16).
Shown below are the graphs of two "solutions" found by shooting. The value of y(1) = y0 differs by only 1e-30 between the two:

It is remarkable how bad it is.
The problem here described can be seen for much higher values of b as mentioned above.

@oceanxinlie What made me think that I was on the right track was the success for this simplified version, where the ode is simply

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

and will be considered on the interval x = 0..1.
Since we only accept y(x)<=0 this abs version should be fine.
Your problem is

ode:=diff(y(x),x,x)=-a*abs(y(x))^0.337+b*x;

on x = 0..2945 with a and b as given by you and boundary conditions y(0)=0, D(y)(2945) = 0.
Your problem can by a change of variables y(x) = A*v(z), x=L*z, be transformed into
 

ode1:=diff(v(z),z,z)=-abs(v(z))^0.337+B*z;

on the interval z = 0..1 and with v(0)=0, D(v)(1)=0, where B = 0.549269696597376e-5.
Thus your problem is close to the my toy problem except for the coefficient B (and the unimportant difference between 0.337 and 1/3). That B is so small seems to be a problem. Right now I don't know why, but will continue looking at it for a while.
##
Here follows the code for the toy problem. It uses a simple shooting method and shoots from x = 1 trying to make y(0) = 0 by varying the value of y(1) = y1. It uses fsolve with a guessed value for y1 found by looking at the animation provided by the code.
dsolve with the parameters option is used.
 

restart;
ode:=diff(y(x),x,x)=-abs(y(x))^(1/3)+x;
plots:-shadebetween(-x^3,0,x=0..1,color=gray); bg:=%:
res:=dsolve({ode,y(1)=y0,D(y)(1)=0},numeric,parameters=[y0],output=listprocedure);
Y:=subs(res,y(x));
Q:=proc(y0,b::truefalse:=true) if not y0::realcons then return 'procname(_passed)' end if;
  res(parameters=[y0]); 
  if b then Y(0) else res end if
end proc;  
plots:-odeplot(Q(-0.1,false),0..1); #Testing Q
plots:-animate(plots:-odeplot,['Q(y0,false)',[x,y(x)],0..1],y0=-0.2..-0.1,background=bg);
fsolve(Q,-0.116); #Finds the value of y(1) for which y(0)=0.
res(parameters=[%]); #Setting the parameter y1.
plots:-odeplot(res,[x,y(x)],0..1,caption="The solution");

@oceanxinlie I'm pretty sure I now have a solution. I'll be gone for a few hours, but will report back then.

Use VectorCalculus:-Jacobian as in
 

VectorCalculus:-Jacobian(L1, L2);
J:=unapply(%,L2);

where L1 is the list of right hand sides and L2 is the list of names (variables).
Follow up with the unapply line and then J will be a function of the variables.

You need pdsolve instead of solve.

@oceanxinlie I didn't use the expression for u that you use above.
I used

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

and remarked at the bottom the note: 
"This 'desperate' attempt to get a positive solution by rewriting the ode ends up just giving the same as the one where u = 0."
So that means I discarded that idea!
Now your expression u is
 

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

but that won't work because (-y(x))^0.337  will give imaginary results for y(x) > 0.
Your previous expression

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

made more sense.
When I asked if y(x) can change sign, I didn't mean to pose a mathematical question, but rater this: Physically (apart from the actual mathematical mode) can y(x) change sign?
It appears that you actually answer 'No' when you say:
"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  "

But (-y(x))^0.337 = abs(y(x))^0.337 for y(x) <=0, so in principle you might as well use
 

ode1:=diff(y(x),x,x)+0.00003019*abs(y(x))^0.337=9.542*10^(-13)*x; 

which, however, will give you a positive solution y(x), which is actually roughly the same as -y(x) from

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

except for the sign. ode2 gives a negative solution as seen in my answer.

If y(x) solves ode1 then z(x) = -y(x) will solve
-diff(z(x),x,x)+0.00003019*abs(z(x))^0.337=9.542*10^(-13)*x;
i.e.
diff(z(x),x,x)-0.00003019*abs(z(x))^0.337=-9.542*10^(-13)*x;

The right hand sides have different signs, but they are small in absolute value.
The maximal relative difference of the solutions of ode1 and ode2 is 1.7e-5.

@oceanxinlie Is it expected that y(x) actually changes sign?
The immediate problem encountered when trying dsolve on your revised ode is that the piecewise expression isn't differentiable at zero:

ode:=diff(y(x),x,x)-0.00003019*piecewise(y(x)>=0,y(x)^0.337,-((-y(x))^0.337))=9.541*10^(-13)*x;
dsolve({ode,bcs},numeric,method=bvp[midrich]); # ERROR

plot(piecewise(y>=0,y^0.337,-(-y)^0.337),y=-1..1);
diff(piecewise(y>=0,y^0.337,-(-y)^0.337),y);

We notice in the numerical solution that R doesn't vary much, between 0.14 and 0.1428.
Thus it might be illuminating to look at the special case where R(theta) in the ode for mu is replaced by a constant R0.
It turns out that this ode with the initial conditions mu(0) = mu0, mu'(0) = 0 has a solution satisfying mu'(Pi) = 0 provided only that R0 and mu0 satisfy the inequality 
-mu0-(1/2)*ln(R0)+(1/2)*ln(2) < 0.
 

restart;
ode0:=R0 = 2*exp(-2*mu(theta))*(1-(diff(mu(theta), theta, theta))-cot(theta)*(diff(mu(theta), theta)));
res:=dsolve(ode0);
resA:=allvalues(res);
## The second result equals the first if _C1 is replaced by -_C1:
resA[1]-eval(resA[2],_C1=-_C1);
## Thus we need only consider the first result resA[1].
###### Graphical experimentation:
eval(rhs(resA[1]),{R0=0.14,_C2= -1,_C1=1});
plot(eval(rhs(resA[1]),{R0=0.14,_C2=-.1,_C1=1}),theta=0..Pi);
plots:-animate(plot,[eval(rhs(resA[1]),{R0=0.14,_C2=-2}),theta=0..Pi],_C1=.9..1.1);
###### Now look at mu'(theta):
diff(rhs(resA[1]),theta);
eval(%,{_C1=1}); #Seems to be a good idea
limit(%,theta=0,right) assuming _C2<0,R0>0;
limit(%%,theta=Pi,left) assuming _C2<0,R0>0;
## So the following expression satisfies mu'(0) = mu'(Pi) = 0 as long as _C2<0 and R0>0:
sol1:=eval(rhs(resA[1]),_C1=1);
## We want mu(0) = mu0 (at least as a limit from the right):
limit(sol1,theta=0,right) assuming _C2<0,R0>0;
c2:=solve(%=mu0,{_C2}); # using that mu(0)=mu0
sol:=eval(sol1,c2); ## Solution for general R0 and mu0
For _C2 to be negative we need R0 > 1/5000:
ineq:=rhs(op(c2))<0;
## The parameters in your case:
params:={R0=0.1428,mu0=ln(100)}; # The "final" value from numerics is taken for R0.
evalf(eval(ineq,params)); #OK satisfied
eval(sol,params); #See the actual solution
plot(eval(sol,params),theta=0..Pi,labels=[theta,mu]); 
plot(eval(diff(sol,theta),params),theta=0..Pitheta=0..Pi,thickness=3,size=[800,default],labels=[theta,diff(mu(theta),theta)]) 

I get no errors in that whole worksheet.
interface(version);
  `Standard Worksheet Interface, Maple 2017.3, Windows 10, September 27 2017 Build ID 1265877`
kernelopts(version);
  `Maple 2017.3, X86 64 WINDOWS, Sep 27 2017, Build ID 1265877`

@Thomas Richard It appears to me that the `*` used and vartheta instead of theta is not a problem affecting the solution at all.

So what in your opinion should C, I, and R be?

There is no attachment!

I don't see any attachment.