@Preben Alsholm

In this part

restart;
Digits:=15:
with(plots);
R0 := 10^(-5);
C := 1;
m := 1/100;
u1 := 1000;
sys := R(x) = 2*(-r(x)*diff(r(x), x)*cos(x)+r(x)^2*sin(x)+diff(r(x), x)^2*sin(x)-sin(x)*r(x)*diff(r(x),x,x))/(r(x)^4*sin(x)),
diff(R(x),x,x)+cot(x)*diff(R(x), x) = -(1/2)*r(x)^2*(R0^2-R(x)^2)-r(x)^2*C^2*m^2*exp(-2*m*A(x))/(2*u1),
diff(A(x), x) = r(x);
## I prefer this version:
SYS:=solve({sys},{diff(r(x),x,x),diff(R(x),x,x),diff(A(x),x)});
## I set the condition on A(x) at x = Pi, but of course we don't know its value.
## We let it be a parameter here called APi.
##
condA := R(Pi) = 1/500^2, r(Pi) = 500, D(R)(Pi) = 0, D(r)(Pi) = 0, A(Pi) = APi;
## Now use the parameters option in dsolve/numeric:
res:=dsolve(SYS union {condA},numeric,parameters=[APi],output=listprocedure);
## Extracting the procedure for finding A(x) for given x:
Aproc:=subs(res,A(x));
## Setting as an example the parameter APi to 10^4:
res(parameters=[10^4]);
## The plots:
plots:-odeplot(res,[x,r(x)],0..Pi);
plots:-odeplot(res,[x,R(x)],0..Pi);
plots:-odeplot(res,[x,A(x)],0..Pi);

you actually find the solution, do you not? The only problem is that we don't know A(Pi) and we can assume it is 10^4, for example. Then we get the plots and solutions for wanted functions, even with a singularity.

I had and issue with it though: I need a plot for original T(x), not A(x), and I get an error: 1.mw

As for the second part, depending on the APi, i. e. A(Pi), there is a singularity at x near to 1. It's a pity that we can't put this singularity outside the interval x=0..Pi, varying A(Pi). Seems kinda strange, as if the system doesn't have a solution without a singularity if the boundary conditions are given at Pi, and yet my professor got it.

As I mentioned before, he used a method of iterations: at first solving the system with T(x)=0, then finding r(x) => T(x) and solving the system with this T(x), repeating 2-3 times. Surprisingly enough, he didn't get a singularity for r(x) from the system without T(x), as you did in

SYS0:=remove(has,eval(SYS,exp=0),A(x));
res0:=dsolve(SYS0 union {R(Pi) = 1/500^2, r(Pi) = 500, D(R)(Pi) = 0, D(r)(Pi) = 0},numeric);
odeplot(res0,[x,r(x)],0..Pi); #Notice the value 0.99622186
odeplot(res0,[x,R(x)],0..Pi);

Although I don't have his worksheet and can be mistaken. However, the final solution for r(x) and T(x) doesn't contain any singularities, that I know.

I've tried to do something similar, but I'm stuck in defining T(x): when I solve the system with a new T(x), Maple doesn't seem to recognise it as a function of x. I think I again have a mistake in the syntax. That is a subject for another question, but maybe you know where I'm mistaken: 2.mw

Thank you for your time and work anyway, you really helped me.