@JoeyJoe Now that you are going to actual data the question is:
Is the ode correct? It appears to me that the dimension of the left hand side of

ode  := diff(1/r(phi),phi$2) + 1/r(phi)  = G*M/h^2 +3*G*M/((c^2)*r(phi));

is 1/L, where L stands for length. The dimension of the right hand side's last term is
(L^3/(M*T^2))*M / ((L/T)^2*L) = 1, i.e dimensionless.
Here T and M stand for time and mass, respectively.
Thus the dimensions of the two sides don't match!
(Edited).

Taking Rouben's reply into account it appears that the first term on the right has dimension 1/L.
Thus a wild guess is that the ode ought to be
 

ode := diff(1/r(phi),phi$2) + 1/r(phi) = G*M/h^2 +3*G*M/((c^2)*r(phi)^2);

With that change I get the value expected 0.33625 in the first case for the differences deltaphi.
Thus I could be right: you forgot a square.
 

@Aaeru Michi The reason we need Aproc(0) = 0 is quite simply because
A(x) = int( r(y), y=0..x).

That's how it all started, remember? Thus A(0) = int( r(y), y=0..0) = 0.
Aproc is just my name for the numerical procedure computing A(x) for given x.
##
You defined T (and I followed you in that) by
T:=exp(-2*m*Aproc(x)); 
Thus T is not a procedure, T is an expression in x. Compare the expression sin(x) with the function (in Maple procedure) sin.
So you cannot apply T to x and expect something sensible. Just as you should not do sin(x)(x).
Notice that I didn't do that. I wrote
plots:-odeplot(res,[x,T],1e-100..0.1,thickness=3);
without an x on T.
##
Finally, you wrote after quoting my comment that "The parameters are not irrelevant":
'I just wanted to solve the system "myself" with random parameters, to "naturally" come to the correct answer, varying the parameters.'

I understand well your wish to solve the system yourself, but is that what you are doing?

@JoeyJoe When I had finished my reply I saw your own reply. Looks good!

@JoeyJoe Here is a numerical solution along the lines of my previous answer, but also answering the perihelia question.
NOTE: I don't find 0.33625, but 0.152650 for the differences.

#############
Extra Note: If the ode below is changed to the dimensionally correct version
 

ode := diff(1/r(phi),phi$2) + 1/r(phi) = G*M/h^2 + 3*G*M/(c^2*r(phi));

then we will find the expected number 0.33625.
This is not changed below, but is easily done.
################

restart;
params:={M = 1,G = 1,h = 1, c=8};
ode  := diff(1/r(phi),phi$2) + 1/r(phi) = G*M/h^2 + 3*G*M/(c^2*r(phi));
## Using output=listprocedure this time:
res := dsolve({eval(ode,params),r(0)= 2/3,D(r)(0) = 0},numeric,output=listprocedure);
## Extracting the procedure that computes r'(phi).
## For minimal distance from the origin we need r'(0) = 0:
R1:=subs(res,diff(r(phi),phi));
## Plotting
plots:-odeplot(res,[r(phi)*cos(phi),r(phi)*sin(phi)],0..10*Pi,scaling=constrained,labels=[x,y],numpoints=1000);
## An animation in phi:
plots:-odeplot(res,[r(phi)*cos(phi),r(phi)*sin(phi)],0..10*Pi,scaling=constrained,labels=[x,y],numpoints=1000,frames=50);
## r as a function of phi:
plots:-odeplot(res,[phi,r(phi)],0..10*Pi);
## We use that graph to help fsolve find the location of the minima:
phi1:=fsolve(R1,5..7);
phi2:=fsolve(R1,11..14);
phi3:=fsolve(R1,18..20);
phi4:=fsolve(R1,25..27);
evalf(phi2-phi1-2*Pi);
evalf(phi3-phi2-2*Pi);
evalf(phi4-phi3-2*Pi);

The animation for the fun of it:

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