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

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