I tried in Windows 10 with Maple 2017.3 64 bit:
 

evalf(Int(x*(1-2*x^(3/10))^(10/3),x=0..1));

That worked fine, but with the exponent 10./3 I got the crash.

@tomleslie As I tried to show in my second answer in the link to the OP's first version of this problem the problem doesn't have any solution defined in the interval 0..Pi.
I wrote:
I think the bad news is that there will necessarily be a singularity in the interior when using conditions at x = Pi.
Thus since int(r(y),y=0..x) appears in the original formulation the conclusion is that the problem has no solution.
In the following I use the revised values you gave in your second worksheet.

Those revised values are same as in this repost: C = 1, u1 = 1000, m=1/100, R0=10^(-5).

@Carl Love I don't see a spam button when viewing this one:
https://www.mapleprimes.com/posts/208776-Avoiding-Scurvy-On-A-Budget

@JoeyJoe Doing just as with the toy parameters and using Rouben's explanation of the meaning of h, you can do:
 

restart;
Digits:=15:
params0:={M = 1.99*10^30,G = 6.67*10^(-11), c=3.0*10^8,r0=4.6*10^10};
#Rouben Rostamian's estimate of h:
h0:=evalf(eval(r0*2*Pi*r0/(88*24*60*60),params0)); #88 days rather than 83
#Using an average orbital velocity of 170496000/3600 m/s found on the internet:
r0*170496000./3600;
hm:=eval(%,params0);
#Using the average of those as a rough estimate:
h1:=(hm+h0)/2;
params:=params0 union {h=h1};
ode  := diff(1/r(phi),phi$2) + 1/r(phi) = G*M/h^2 + 3*G*M/(c^2*r(phi)^2);
## Using output=listprocedure:
res := dsolve(eval({ode,r(0)= r0,D(r)(0) = 0},params),numeric,output=listprocedure,abserr=1e-15,relerr=1e-13,maxfun=10^6);
## 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);
#r as a function of phi:
plots:-odeplot(res,[phi,r(phi)],0..10*Pi,size=[1200,default]);
## We use that graph to help fsolve find the location of the minima:
phi1:=fsolve(R1,3..5);
phi2:=fsolve(R1,8..10);
phi3:=fsolve(R1,15..17);
phi4:=fsolve(R1,21..23);
evalf(phi2-phi1-2*Pi);
evalf(phi3-phi2-2*Pi);
evalf(phi4-phi3-2*Pi);
phi1:=fsolve(R1,3..5);
### ## Attempting a more automatic approach:
for i from 1 to 15 do 
  phi2:=fsolve(R1,phi1+2*Pi);
  delphi:=evalf(phi2-phi1-2*Pi);
  print(delphi);
  phi1:=phi2
end do:

I get 9.569995*10^(-7) for delphi, but with another and more accurate h you will find something better I hope.

Actually I think I already answered that question in this
https://www.mapleprimes.com/questions/223373-How-Do-I-Put-An-Integral-With-Variable

So why ask again?

@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]?