@Aaeru Michi Using the parameters in figure 1 I can reproduce the graphs using the parameters option in dsolve as given in the code in my second answer. It is clearly not irrelevant what the parameters are.
For clarity I bring the whole code. I have simplified Q since the purpose no more is to show that a singularity exists in the interior, but now we need to find a value for APi that makes Aproc(0) = 0, or as it turns out Aproc(1e-150) = 0 which should be OK.

restart;
Digits:=15:
with(plots):
R0 := 10^(-5);
C := 1;
m := 0.15; ## As in figure 1
u1 := 100; ## As in figure 1
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);
SYS:=solve({sys},{diff(r(x),x,x),diff(R(x),x,x),diff(A(x),x)});
condA := R(Pi) = 2/447^2, r(Pi) = 447, D(R)(Pi) = 0, D(r)(Pi) = 0, A(Pi) = APi; #As in figure 1
res:=dsolve(SYS union {condA},numeric,parameters=[APi],output=listprocedure,abserr=1e-13,relerr=1e-11,maxfun=10^6); #Added stricter tolerances than default.
Aproc:=subs(res,A(x));
## Preliminary testing:
res(parameters=[2000]);
plots:-odeplot(res,[x,r(x)],0..Pi);
plots:-odeplot(res,[x,R(x)],0..Pi);
plots:-odeplot(res,[x,A(x)],0..Pi);
## The simplified Q procedure:
Q:=proc(aa) local out,xx;
   xx:=1e-150; # A substitute for zero.
   if not type(aa,realcons) then return 'procname(_passed)' end if;
   res(parameters=[aa]); 
   try 
     out:=Aproc(xx);
   catch:
     out:=-10^3
   end try;
   out
end proc:
## Testing Q and plotting:
Q(2000);
plot(Q(a),a=1300..2000,adaptive=false,numpoints=10,style=point,symbolsize=20);
plot(Q(a),a=1390..1450,adaptive=false,numpoints=15,style=point,symbolsize=20,view=-100..50);
plot(Q(a),a=1400..1410,adaptive=false,numpoints=15,style=point,symbolsize=20);
## We see that we must have a zero for Q between 1400 and 1403.
## We use Student:-NumericalAnalysis:-Secant instead of fsolve:
APival:=Student:-NumericalAnalysis:-Secant(Q(a,1e-100),a=[1400,1403]);
res(parameters=[APival]); #Setting APi to APival
## The final plots:
plots:-odeplot(res,[x,r(x)],0..Pi);
plots:-odeplot(res,[x,R(x)],1e-100..Pi,thickness=3);
plots:-odeplot(res,[x,R(x)],0..0.1,thickness=3);
T:=exp(-2*m*Aproc(x));
plots:-odeplot(res,[x,T],1e-100..0.1,thickness=3);

Here are two of the plots (r and T):

@Aaeru Michi I would like to see his solution.

@Aaeru Michi Yes, dsolve/numeric/bvp uses iteration just starting from a very rough first guess unless you have some prior knowledge of the solution you can hand to it. It is just like using Newton's method. (In fact it probably uses that method at least in part of the process).
I suggest that you look at my second answer below, which concludes that the whole problem has no solution. A bold statement, which takes into account that you need the solution to exist on the whole interval 0..Pi.
The problem is in fact not a singularity at the endpoints, but a singularity in the interior.
It makes me think of a very simple ode which exhibits a singularity very nicely:
 

restart;
ode:=diff(y(x),x)=y(x)^2;
dsolve({ode,y(0)=1}); # singularity at x = 1
## Now numerically and first using dsolve/numeric/bvp on 0..Pi:
resBVP:=dsolve({ode,y(0)=1},numeric,method=bvp,range=0..Pi);
## Using the fact that it is an initial value problem:
resINI:=dsolve({ode,y(0)=1},numeric);
plots:-odeplot(resINI,[x,y(x)],0..Pi);

It appears that you have an empty file, 0 kb. Thus nothing can be recovered.

After looking once more at your image above I wondered why there are two different notations for derivatives w.r.t. theta, one being a prime (as in r') the other being written as a partial derivative of R w.r.t. theta.
In your worksheets you treat them the same, i.e. ordinary derivatives w.r.t. theta (OK you use x instead of theta).
Do you have a link to the source for the image, which may throw some light on this?
Maybe the meaning is that the first line defines the expression R and in the second line the partial of R w.r.t. theta is supposed to differentiate the expression R as a function only of the explicitly appearing theta, i.e. that appearing in sin(theta) and cos(theta)?

I just noticed that the image has r[b] in the integral, not r. You have r in your worksheet. What is r[b]?

@Aaeru Michi By default dsolve/numeric/bvp finds an approximate solution itself. That one is based entirely on the boundary conditions. In your case with

R(Pi) = 1/500^2, r(Pi) = 500, D(R)(Pi) = 0, D(r)(Pi) = 0, A(0) = 0;

it uses the approximate solution

approxsoln=[A(x) = 0, R(x) = 1/250000, r(x) = 500]

As you will notice this approximate solution just consists of constant functions, but they certainly satisfy the boundary conditions. Takíng the ode system itself into account you may be able to do better by explicitly giving an option of the form

approxsoln = [A(x)= Ac(x), R(x)=Rc(x), r(x)=rc(x)]

where each of Ac, Rc, and rc are concrete functions of x that you think might get the numerical solution started well.
Although I have tried I haven't been successful with my guesses for Ac, Rc, or rc.

I tried solving for r(x) using the ode with diff(r(x),x,x) and just setting R(x) = 1/500^2. Specifically I did this:
 

ode_r:=diff(r(x), x, x) = -(1/2)*r(x)^3*R(x)+r(x)-diff(r(x), x)*cos(x)/sin(x)+diff(r(x), x)^2/r(x);
res:=dsolve({eval(ode_r,R(x)=1/500^2),r(Pi)=500,D(r)(Pi)=0},numeric,output=listprocedure);
res(0); #Problem very, very close to 0.
rp:=subs(res,r(x));
evalf(Int(rp,0..Pi));
plot(Int(rp,0..x),x=0..Pi,numpoints=10,adaptive=false);
plots:-odeplot(res,[x,r(x)],0..Pi);

That was done to get an idea of a guess for r(x) and A(x).  I haven't had any successes so far, though.

@Markiyan Hirnyk You are quite right. But if eager0626 wants to see the roundoff then plotting the difference makes sense.

@torabi You have a system which is a variant of the famous Lorenz equations. But you ought to explain for what values of t you measure y(t) for different values of A.
In the meantime here is a nice picture to look at. A is 25.

It is generated by the following using your code except for changing A:=1 to A:=25.
 

plots:-odeplot(dsol5,[x(t),y(t),z(t)],10..110,numpoints=10000,caption="A = 25");

@Christopher2222 

You can try this as an example:
 

restart;
type(a[1],symbol);
type(a[1],name);

@John Fredsted Another simple way using that exp(1/i) > 1 for all i >=1 (thus P is an increasing function of n):

P := product(exp(1/i),i = 1..n);
fsolve(P=100,n);

The output from fsolve is 55.64520627, so the answer is 56.

@_Maxim_ Just a comment about the last point. Notice the difference between these two calls, where the sequence defined as a is evaluated by evalf at the setting of Digits, and the second where evalf is given the explicit sequence Pi,7 whereby Pi is evaluated at Digits=7:

restart;
Digits:=10:
a:=Pi,7;
evalf(a);
evalf(Pi,7);

@tomleslie There is no doubt that the system dsys3 satisfies the requirements for existence and uniqueness.
Solving for the second derivatives gives right hand sides that are just a sum of terms including besides sin(t) and constants only polynomials of order 3 in the dependent variables and their derivatives. 

But how large the maximal interval of definition of the solution is, is certainly not clear. It may be (1) finite, and (2) very small. Thus it might not even include  t = 1.

I notice that your equations eq3 and eq4 contains integrals. This one is from eq4:
int(cosh(4.6941*x)-cos(4.6941*x)-.9825251551*sinh(4.6941*x)+.9825251551*sin(4.6941*x), x = 0 .. .999)

and it just evaluates to .8369623967. The similar integral found in eq3 just evaluates to .2866222667.

Is this as intended?

@Carl Love Since the values in the output = array(...) are the values of the independent variable, it doesn't make much sense to ask for t=1 three times.