Although I tried many times before I wrote the question above and had no success, I finally had success after taking a short break.
This is being written using Firefox.

In previous days I had run into the same problem with the new Firefox, but trying once or twice more got me in.
I'm still wondering what is going on. I have no problems with other websites that I visit.

This could happen e.g. in the following circumstance:
 

seq(diff(f(x),x$i),i=1..3); #OK
seq(diff(f(x),x$i),i=0..3); # Your error
seq(diff(f(x),[x$i]),i=0..3); # OK with list
diff(f(x),[]); # OK just f(x)

@mrmusic1994 Besides what Kitonum has already replied, notice that if k < 0 then the "terminal velocity" is infinity.
Using the notation from my answer we would find that Vfinal is infinity. Thus no negative k can solve you problem!

Vfinal:=limit(V,t=infinity) assuming k<0;

It is equally easy to see that k = 0 won't solve your problem either. Do a separate calculation:
 

de0 := 100*diff(v(t),t)=100*9.81;     
dsolve({de0,v(0)=0});

Again the "terminal velocity" is infinity.

@edmundac I tried your worksheet in Maple 2017.3 as well as in Maple 17.02 on my Windows 10 64-bit computer.

I ran the code once. Then I ran the same worksheet once more from the very beginning, i.e. starting with restart. I didn't find any problem. So I wonder what you did?
Did you do something differently the second time?

@killercrew You can change the outer loop to
 

for i from 0 to N-1 do

and still keep the inner as I wrote. The change in the outer loop just means a renumbering, so that if X and Y are defined as Arrays then without the outer loop change they would be
 

X:=Array(1..M,0..N,datatype=float);
Y:=Array(1..M,0..N,datatype=float);

whereas with the outer loop changed they would be
 

X:=Array(0..M-1,0..N,datatype=float);
Y:=Array(0..M-1,0..N,datatype=float);

The datatype option I added. That has nothing to do with the numbering.

@vv Since y' = sqrt(2-2*ln(y^2))*y  it follows that as long as y > 0 we have y' >=0. Thus y(x) is nondecreasing for all x for which the solution is defined and y(x) > 0. The ode has another problem at y = 0, but that is of no concern here.

But I agree that there is no bug. No claim is made by the help page that a maximal solution is given.
 

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: