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.

@vv Just a comment on infolevel: With  infolevel[Basis]:=3:  (or 2) we do get information, but it may not be of interest to codell.
Since the procedure remembers the information will only come the first time it is run with the given input.

I have never heard of acximetric tolerance. Is it spelled correctly?
You could be a little more specific and provide an actual example.
When you say " For example +10/-5" I still don't get what you are talking about.

I read at the top of the pdf file containing the problems as formulated by your teacher:
"Studentene får hver sin oppgave, som skal løses selvstendig."'

Does that not exclude getting help from MaplePrimes? Ask your teacher.

@tomleslie I have noticed before that these problems come up too often for students in Denmark.
I cannot read the uploaded file either.

Personally, I haven't had these problems in years although my computer is made for use with Danish characters and bought in Denmark.
Years ago when I occasionally had problems, I started a habit of never (well, almost never) saving a worksheep with output.
I routinely remove all output from the worksheet before saving it by going to
Edit/Remove Output/From Worksheet.
Another thing is that I never use 2D math input or fancy text. My worksheets are rather plain, but I do use the Standard Interface (with Maple Input and Worksheet mode).
##
acer has been able to rescue some people in the past.

@Markiyan Hirnyk Where do you find HINT=sum documented?
HINT=xxx works the same:
 

restart;
PDE := (y-u(x, y))*diff(u(x, y), x)+(u(x, y)-x)*diff(u(x, y), y) = x-y;
infolevel[pdsolve] := 3: 
ans := pdsolve(PDE, HINT = xxx);