@tomleslie A few hours ago I filed a bug report about the following example from the help page for LPSolve:
 

restart;
Optimization:-LPSolve(2*x+5*y, {3*x-y=1, x-y<=5}, assume={nonnegative, integer});

It produces a crash (loss of kernel connection) on my computer in Maple 18.02, 2015.2, 2016.2, and 2017.3, but NOT in Maple 17.02 where it gives a correct result: [12, [x = 1, y = 2]]

Do you get a the crash on your computer in the Maple versions you have?

@tomleslie For Maple 2015.2 I get

restart;
kernelopts(version);
Optimization:-LPSolve( 3*x__1+14*x__2+18*x__3+6*x__4+2*x__5,
                       {3*x__1+5*x__2+6*x__3+2*x__4+x__5 <= 10},
                       binaryvariables=[x__1, x__2, x__3, x__4, x__5],
                       maximize = true
                      );
  Maple 2015.2, X86 64 WINDOWS, Dec 20 2015, Build ID 1097895
Warning, problem appears to be unbounded
    [0, [x__1 = 0, x__2 = 0, x__3 = 0, x__4 = 0, x__5 = 0]]

#########
Maple 18.02:

restart;
kernelopts(version);
Optimization:-LPSolve( 3*x__1+14*x__2+18*x__3+6*x__4+2*x__5,
                       {3*x__1+5*x__2+6*x__3+2*x__4+x__5 <= 10},
                       binaryvariables=[x__1, x__2, x__3, x__4, x__5],
                       maximize = true
                      );
   Maple 18.02, X86 64 WINDOWS, Oct 20 2014, Build ID 991181
Warning, problem appears to be unbounded
    [0, [x__1 = 0, x__2 = 0, x__3 = 0, x__4 = 0, x__5 = 0]]

#############
But in Maple 17.02:

restart;
kernelopts(version);
Optimization:-LPSolve( 3*x__1+14*x__2+18*x__3+6*x__4+2*x__5,
                       {3*x__1+5*x__2+6*x__3+2*x__4+x__5 <= 10},
                       binaryvariables=[x__1, x__2, x__3, x__4, x__5],
                       maximize = true
                      );
    Maple 17.02, X86 64 WINDOWS, Sep 5 2013, Build ID 872941
    [26, [x__1 = 0, x__2 = 0, x__3 = 1, x__4 = 1, x__5 = 1]]

##########################################
NOTE ADDED:
I thought that the build ID number would be the same for the same final release for (say) Windows 64 bit. By "final" I exclude beta versions.

@tomleslie I tried your code in Maple 2017.3 and in Maple 2016.2 and got the warning about unboundedness and the result
[0, [x__1 = 0, x__2 = 0, x__3 = 0, x__4 = 0, x__5 = 0]]

For Maple 2017.3:

kernelopts(version);
  Maple 2017.3, X86 64 WINDOWS, Sep 27 2017, Build ID 1265877

For Maple 2016.2:

kernelopts(version);
  Maple 2016.2, X86 64 WINDOWS, Oct 21 2016, Build ID 1174570

###############

But in Maple 17.02 I got the result you report.
The same goes for Maple 12, 15, and 16.
################
In Maple 18 and Maple 2015 it is the same as in Maple 2017.3.

The problem as formulated can be solved by NLPSolve.
 

f:=3*x__1+14*x__2+18*x__3+6*x__4+2*x__5;
cstr:={3*x__1+5*x__2+6*x__3+2*x__4+x__5 <= 10};
res:=Optimization:-NLPSolve(f,cstr, x__1 = 0 .. 1, x__2 = 0 .. 1, x__3 = 0 .. 1, x__4 = 0 .. 1, x__5 = 0 .. 1, maximize = true);
eval([f,cstr],res[2]); #Check

This doesn't solve the binary problem, but neither would the command you are using (if it had worked).

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?