I think you should use dsolve,numeric.
I have chosen to use the parameters option, but you don't have to do that.
 

restart;
ode:=diff(u(t),t$2)-0.1*(1-64*u(t)^2)*diff(u(t),t$1)+16*u(t)=0;
Q1:=v/4*sin(4*t)+epsilon*(1/8*(v-v^3)*t*sin(4*t)+v^3/128*(cos(4*t)-cos(12*t)));
##
res:=dsolve({ode,u(0)=0, D(u)(0)=v},numeric,parameters=[v]);
##
res(parameters=[2]); #Setting the parameter v to 2 in res
A:=plots:-odeplot(res,[t,u(t)],0..10,legend="Numerical solution");
B:=plot(eval(Q1,{v=2,epsilon=0.1}), t=0..10, u=-0.5..0.5,color=blue, legend="Regular Perturbation Expansion");
##
plots:-display([A,B],title="Comparing the regular perturbation expansion to the numerical solution");

## For the fun of it we can animate in epsilon with v fixed or vice versa:
 

Q:=proc(v1,eps,{range::range:=0..10}) local A,B;
  if not [v1,eps]::list(realcons) then return 'procname(_passed)' end if;
  res(parameters=[v1]);
  A:=plots:-odeplot(res,[t,u(t)],range,legend="Numerical solution");
  B:=plot(eval(Q1,{v=v1,epsilon=eps}), t=range,color=blue, legend="Regular Perturbation Expansion");
  plots:-display(A,B,title="Comparing the regular perturbation expansion to the numerical solution")
end proc:
  
Q(2,0.1);
Q(2,0.1,range=0..2);
plots:-animate(curry(Q,2),[epsilon,range=0..1],epsilon=0..1); # v fixed at 2
plots:-animate(rcurry(Q,0.1),[v,range=0..5],v=1..2); # epsilon fixed at 0.1

Never use r and r[aaa] in the same command; here aaa can be "anything".
You can use r[i], r[bc], r[tc] in the same command though.
The background is that as soon as the assignment r[aaa]:=qqq is made the name r now refers to a table.
Try

restart;
r[23]:=17;
eval(r);

You will see that r refers to the table table([23 = 17]), which means that the indexed name r[23] has the value 17.
Although you yourself didn't make any assignments, this fact can create havoc anyway. It isn't always clear why.
The use of double underscores, as in r__bc doesn't have any similar problem r__bc prints just as r[bc] does.
 

int(r*r__bc*r__tc, r = r__bc .. r__tc);

I tried this:
 

restart;
#k=2:
omega:=u/h:
psi:=v/h:
t:=(add(a[j]*x^j,j=0..2)+a[3]*sin(omega*x)+a[4]*cos(omega*x)+a[5]*sin(psi*x)+a[6]*cos(psi*x)):
F:=diff(t,x):
G:=diff(t,x,x):
p1:=simplify(eval(t,x=q+h))=y[n+1]:
p2:=simplify(eval(F,x=q))=f[n]:
p3:=simplify(eval(F,x=q+h))=f[n+1]:
p4:=simplify(eval(F,x=q+2*h))=f[n+2]:
p5:=simplify(eval(G,x=q))=g[n]:
p6:=simplify(eval(G,x=q+h))=g[n+1]:
p7:=simplify(eval(G,x=q+2*h))=g[n+2]:
vars:= seq(a[i],i=0..6):
Cc:=eval(<vars>, solve({p||(1..7)}, {vars})):
for i from 1 to 7 do
	a[i-1]:=Cc[i]:
end do:
Cf:=t:

K:= collect(combine(simplify(eval(Cf,x=q+2*h),size),trig),{y[n+1],f[n],f[n+1],f[n+2],g[n],g[n+1],g[n+2]},factor):

## Notice assignments and for gamma name changes:
alpha[1]:=simplify(coeff(K,y[n+1]));
beta[0]:=simplify(coeff(K,f[n]),size);
beta[1]:=simplify(coeff(K,f[n+1]),size):
beta[2]:=simplify(coeff(K,f[n+2]),size):
gamma0:=simplify(coeff(K,g[n]),size):
gamma1:=simplify(coeff(K,g[n+1]),size):
gamma2:=simplify(coeff(K,g[n+2]),size):
length~([alpha[1],beta[0],beta[1],beta[2],gamma0,gamma1,gamma2]);
length~(simplify~(expand~([alpha[1],beta[0],beta[1],beta[2],gamma0,gamma1,gamma2]),size));

The lengths reported in Maple 2017.3 (and in Maple 2016.2) were
[1, 1945, 1765, 2129, 1941, 2039, 2073] and [1, 1381, 1349, 1429, 1377, 1425, 1441].
Thus some reduction in size.

I converted your expression to 1D math input (aka Maple input).
Called it EXPR.
Then I did

length(EXPR);
EXPR1:=simplify(EXPR);
length(EXPR1);

The first length was 15601, the second 3773. So a considerable reduction in size.
I used Maple 2017.3.
MaplePrimes18-02-19simplify.mw

The reason I converted to 1D input was that your worksheet didn't seem to respond to ENTER. Nothing seemed to happen. I just right now tried copying the whole expression and pasting in a fresh worksheet.
Then I could do as described above without converting to 1D input.
###  You had a colon at the end, not a very good idea. You won't see anything!

I looked at the help page for copy and found the following two statements (in Maple 2017):

"If a is not a table, rtable, or record, and a deep copy is not requested, a is returned."
"To obtain a full recursive copy where all contents that are also mutable structures are copied, or to obtain a copy of a module, use the deep option."

Therefore I tried replacing copy(sol) by copy(sol,deep) in both places.
It seemed to do what you want.

## Note: You have Maple 2015. The option works there, but is not documented.
It is as of Maple 2016.
Try in Maple 2015:

showstat(copy);

You will see that the first line is
if 1 < nargs and args[2] = (':-deep') then ...

Maple procedures evaluate to their names.
Thus your procedures named T1 and T2  just evaluate to those names.
But first of all you need to define T1 and T2 properly.
You must use unapply:

T1:=unapply(Q1,m);
T2:=unapply(Q2,m);

The reason is that in doing T1:= m-> Q1, the m in Q1 is not seen at the time of definition, only after a call.
A simple example:

QQ:=x^2;
f:=x->QQ;
##Try
f(5);

It doesn't work. Use unapply or use the direct definition f:= x -> x^2.
 

Use the evalc command.
evalc will assume that all variables are real and will then write the input expression on the form a+I*b.
Simple example:

evalc( exp(I*t) );

Since you indeed can solve successfully for u[3] = 0 you can try continuation.
Change the definition of u[3] to include a factor e.g. called cont and then use the option continuation=cont.
I had success with the following:
 

## All the previous code and then the following instead of the original u[3]:
u[3] := min(max(0, e), 1)*cont;
## Now the rest of your code until dsolve, which is replaced by:
p1:=dsolve({sys}, numeric, abserr = 0.1e-4,continuation=cont,maxmesh=2048);

Your plot:

I don't get an error message from your code. I did, however, replace the print statement with an assignment of the successful loop parameters to an array since the printing takes up an enormous amount of space. There were 854590 successful cases out of a total of 1078110. (!!)

restart;
som:=0;
tot:=0;
A:=Array():  #Added
N:=10;
for a from 4000 to 5000 by 100 do
for b from 100 to 900 by 10 do
for c from 3000 to 3000 by 100 do 
for d from 3000 to 3100 by 10 do
for e from 10000 to 10000 by 100 do
for f from 600 to 700 by 10 do
for p from 1 to 10 by 1 do
tot:=tot+1;
C2:=c*d-N^2*(f^2-b*d);
C3:=a*d-e*f;
C4:=c*e-N^2*(a*f-e*b);
if (C2>0 and C3>0 and C4>0) then
   som:=som+1;
  A(som):=[a,b,c,d,e,f,p]
fi;

od;od;od;od;od;od;od;
##
som;
tot;
A[som]; #Last entry
A[1..5]; # First 5 entries

This runs rather fast.
If in another case you run into exceptional cases resulting in error and just want to skip those you can do as in this rather trivial example:
 

res:=0:
for i from 1 to 10 do
  try
    res:=res+1/(i-3)^2
  catch:
    printf(cat(StringTools:-FormatMessage(lastexception[2..-1]),"\n"));
    next
  end try
end do;
###
res;

You can leave out the printing of the error message. You can also replace catch: with catch "numeric exception":
as in this code:
 

res:=0:
for i from 1 to 10 do
  try
    res:=res+1/(i-3)^2
  catch "numeric exception":
    next
  end try
end do;
##
res;

By (-8)^(1/3) is meant the principal cube root of -1.
The polynomial x^3 + 8  has three roots:  -2, and the pair 1-I*sqrt(3) and 1+I*sqrt(3) .
The principal root is the one with the least (principal) argument in absolute value. In case of equality of absolute values (as here) the one that is positive is chosen.
The principal arguments of the 3 roots are Pi, -Pi/3, and Pi/3 in the order used above.
So the root having argument Pi/3 must be chosen, thus it is 1+I*sqrt(3).

restart;
beta:=simplify( (-8)^(1/3) );
solve(x^3+8=0,x);
argument~([%]);
convert((-8)^(1/3), RootOf);
evalc(%);
simplify(%);

Is this an example of what you have in mind?
 

A:=Array(1..3,1..2,(i,j)->[seq(k,k=i..i^2+j^2)]);
numelems~(A);
numelems(A[2,2]);

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);