The floats in the output are hardware floats.
The very interesting observation by Tom Leslie is illustrated in the following code:
 

restart;
Digits:=10: #Default
X := [seq(0.1e-1*i, i = 0 .. 12)]; 
Y := [0, .36, .56, .67, .73, .76, .78, .79, .79, .80, .80, .80, .80]; 
f := Statistics:-Fit(a-b*exp(c*x+d), X, Y, x, initialvalues = [a = .8, b = 1, c = -50, d = -.3]);
indets(f,hfloat);
fs:=evalindets(f,hfloat,convert,sfloat); # converting to sfloat
indets(fs,hfloat); #none
evalf[8](f);
evalf[11](f);
evalf[11](fs);
## Now constructing a hardware float from a software float:
restart;
Digits:=12: #Changed for illustration only. Digits=10 gives the same.
f:=1.234567890987654321;
fh:=HFloat(f);
type~([f,fh],hfloat);
for i from 3 to 21 do
  forget(evalf);
  i,evalf[i](fh);
end do;
for i from 3 to 21 do
  forget(evalf);
  i,evalf[i](f);
end do;

We see that with the software float f evalf behaves as expected whereas on the hardware float fh evalf has its own (peculiar) behavior.

Edited: You must be using an earlier Maple version than Maple 18.
My first response was based on a check in Maple 2017.
First response:
Your title, ""simplify seems to have a bug here?" couldn't only be based on you not being able to see why this is correct, or is it?
If you suspect a bug, then give us a reason better than that, e.g. an N for which this is false.
##
It should be noticed that the output from simplify seen in your question omits the summation sign and thereby makes the result look nonsensical.
It is OK in your worksheet.  
## Whoops, actually it isn't OK in your worksheet.
## Maple 15, 16, and 17 produce the nonsensical output seen in your question.
Nonsensical because the output refers to the summation index.
Maple 18, Maple 2015, 2016, and 2017 give a sensical answer and an answer that I have no reason to believe is wrong.
## A Maple 2017 worksheet:

 

## Note: the result for res3 holds for NF >= 5, not only NF>5.
 

Download MaplePrimes18-02-10_sum_floor.mw

As stated in the Optimization help page:
"The Optimization package is a collection of commands for numerically solving optimization problems, ... "

thus you cannot use Minimize from that package. But there is no need to either.
(1) Solve cnsts for x.
(2) Replace x in obj with the expression from (1). Call the result f.
(3) f is a polynomial of degree 2 in y, f=A*y^2 + B*y + C. Find the coefficients A, B, and C.
(4) Use the well-known criterion for a minimum of f and its location when there is one.

solve does this with no problem:
 

solve(sin(sqrt(x)),x,allsolutions);

The variable _Z1 stands for an arbitrary integer:
 

about(_Z1);

The short answer is no! The signs of b(x) and a(x) are the ones that are important, not x.

If you look at the help page for arctan you will find this statement:
"For real arguments x, y, the two-argument function arctan(y, x), computes the principal value of the argument of the complex number  x + I y, so  -Pi < arctan(y, x) and arctan(y, x) <= Pi."
 You could try this:
 

arctan(y,x) assuming x>0, y::real;
arctan(y,x) assuming x<0, y>=0;  
arctan(y,x) assuming x<0, y<0;
arctan(y,0) assuming y>0;
arctan(y,0) assuming y<0;
arctan(0,0);

With the default Digits=10 I actually get an error:

Error, (in Optimization:-LPSolve) invalid input parameter to external routine

Raise Digits to 25 and you will get the result you expect. So do
Digits:=25;
before using Maximize.

In doing
evalf( exp(1) );
you are getting a decimal number approximation to exp(1). With Digits at 10 ( the default) that number is 2.718281828 .
You may try
evalf[50](exp(1));
and you get 2.7182818284590452353602874713526624977572470937000 .
Try replacing 50 by e.g. 500.
Maple will not give you an approximation of a number like exp(1) or Pi unless you ask for it either directly as above or indirectly as in
exp(1.);
exp(1.0); #same thing
where I have replaced the integer 1 with the decimal number 1.0 (or just 1. ).
I^2 evaluates to -1 exactly.
But try
I^2.0;
sqrt(2)^2 evaluates to 2 exactly.
But try
sqrt(2.0)^2;
Don't use decimal numbers if you are not satisfied with an approximation.

If you don't need the various statistical information that can be supplied by LinearFit you can use LSSolve:
 

restart;
pts:=[[0,0],[1,13],[2,21],[6,45],[12,54],[15,77]];
X:=evalindets[2](pts,list(realcons),1,op);
Y:=evalindets[2](pts,list(realcons),2,op);
f:=x->a*(x-2)+21; # y-21 = a*(x-2)
d:=f~(X)-Y;
Optimization:-LSSolve(d);

The successful methods are ftoc and ftocms (ftoc stands for The Fundamental Theorem of Calculus).
If you find an antiderivative:

F:=int(sin(K*x)*sin(L*x)*cos(M*x), x);

you will see that
F := (1/4)*sin((K-L-M)*x)/(K-L-M)+(1/4)*sin((K-L+M)*x)/(K-L+M)-(1/4)*sin((K+L-M)*x)/(K+L-M)-(1/4)*sin((K+L+M)*x)/(K+L+M)

thus this is only valid if K-L-M, K-L+M, K+L-M, and K+L+M are all different from zero.
In your case K-L-M = 14-2-12 = 0, so the generic F cannot be used.
### WORKAROUND:
You can use that if a parameter is present in the limits then the option AllSolutions will give all cases.
 

restart;
A:=Int(sin(K*x)*sin(L*x)*cos(M*x), x = 0 .. Pi);
## We may assume that K<>0 since otherwise A = 0 no matter the values of L and M.
B:=IntegrationTools:-Change(A,x=t/K,t); # Assuming that K<>0
f,r:=op(op(2,B));
int(f,r,AllSolutions)/K assuming integer;
C:=simplify(%);
eval(C,{K=14,L=2,M=12});

Your system may have no solution or at least such a solution is very difficult to find for the following reason.
Solving for the highest derivatives gives you a system whose right hand sides has a denominator that will take the value 0 in the open interval (0, L). That just follows from the boundary conditions as you will see here:

restart;
Eq1 := diff(F(eta), eta, eta, eta)-phi1*((diff(F(eta), eta))^2-(diff(F(eta), eta)+diff(G(eta), eta))*(diff(F(eta), eta, eta)))+alpha*(1-phi)^2.5*((diff(F(eta), eta, eta))^2+2*(diff(F(eta), eta))*(diff(F(eta), eta, eta))-(diff(F(eta), eta)+diff(G(eta), eta))*(diff(F(eta), eta, eta, eta, eta))+A*(2*(diff(F(eta), eta, eta, eta))+(1/2)*eta*(diff(F(eta), eta, eta, eta, eta))))-A*phi1*(diff(F(eta), eta)+(1/2)*eta*(diff(F(eta), eta, eta))) = 0; Eq2 := diff(G(eta), eta, eta, eta)-phi1*((diff(G(eta), eta))^2-(diff(F(eta), eta)+diff(G(eta), eta))*(diff(G(eta), eta, eta)))+alpha*(1-phi)^2.5*((diff(G(eta), eta, eta))^2+2*(diff(G(eta), eta))*(diff(G(eta), eta, eta))-(diff(F(eta), eta)+diff(G(eta), eta))*(diff(G(eta), eta, eta, eta, eta))+A*(2*(diff(G(eta), eta, eta, eta))+(1/2)*eta*(diff(G(eta), eta, eta, eta, eta))))-A*phi1*(diff(F(eta), eta)+(1/2)*eta*(diff(F(eta), eta, eta))) = 0;
###
bcs:=F(0)=0,D(F)(0)=1,G(0)=0,D(G)(0)=p, D(F)(L)=0, D(G)(L)=0, F(L)=0, G(L)=0;
###
phi1 := 1.2; alpha := 2; A := 1.5; phi:=0.1; p:=0.1; L:=8;
sys:=solve({Eq1,Eq2},diff~({F,G}(eta),eta$4));
(denom@rhs)~(sys); # The same denominator in both equations
DN:=op(%);
eval[recurse](convert(DN,D),{eta=0,bcs}); # Value 1.690553637*10^9 > 0
eval[recurse](convert(DN,D),{eta=L,bcs}); # Value -9.221201656*10^9 < 0

Even if you change the late coming two boundary conditions won't alter the result of the last two lines: Try using

bcs6:=F(0)=0,D(F)(0)=1,G(0)=0,D(G)(0)=p, D(F)(L)=0,D(G)(L)=0;

instead of bcs in the last two lines above.

Try this:
 

restart;
f:=x->sqrt(2.0)*x;
y := sqrt(3.0); #Digits = 10
for n from 10 to 3 by -1 do
  [evalf[n](evalf[n](sqrt(2.0))*evalf[n](y)),evalf[n](f(y))]
end do;

The two numbers for all n = 3..10 are exactly the same.
And by the way, so are these:

for n from 10 to 3 by -1 do
  [evalf[n](sqrt(2.0)*y),evalf[n](f(y))]
end do;

Go to Tools/Options/Display and uncheck Show equation labels.
Apply Globally or Apply to Session.

restart;
f:=x->x^2-3;
x:=1;
for i from 2 to 5 by 1 do
    x:=evalf(x-f(x)/D(f)(x));
    print(x);
end do;
evalf(sqrt(3));

Besides what Kitonum has already pointed out, you must leave 'continuation' out. There is no continuation parameter in your system.

Furthermore, it appears that the space between Array and the list of numbers is being interpreted as multiplication (as it will be in 2D math input):
It should be:
 

output = Array([0., 0.5e-1, .10, .15, .20, .25, .30, .35, .40, .45, .50, .55, .60, .65, .70, .75, .80, .85, .90, .95, 1.00])

Finally, for plotting you only need to do

with(plots); 
odeplot(r1, [x, f(x)]);

No need for display (you missed an assignment there too, incidentally).
 

The problem is that when the procedure is called as in JD(2018, 1, 11) , then the first that happens is that all occurrences of YY, MM, and DD in the procedure body are replaced with the numbers 2018, 1, and 11 respectively.
The assignments YY := YY-1; MM := MM+12 thereby becomes
2018:=2018-1; 1:=1+12;
But needless to say, you cannot assign to numbers like 2018 or 1.
So change to
 

JD := proc (Y, M, DD)
  local A, B, YY, MM;
  YY:=Y; 
  MM:=M;
  if MM = 1 or MM = 2 then
    YY := YY-1; 
    MM := MM+12
  end if;
  A := floor((1/100)*YY);
  B := 2-A+floor((1/4)*A);
  evalf(floor(365.25*YY)+floor(30.6001*(MM+1))+DD+1720994.5) 
end proc;

Whether your procedure works as intended is another matter.