As we all agree this is just a "cleaning up" business. There is nothing wrong with the solution.
I just tried to make a procedure CleanUp. It assumes that the input is a solution to a linear ode and comes from dsolve.
It hasn't been tested thoroughly.
For what it is worth here it is:
 

CleanUp:=proc(sol::equation) local t,cs,Hom,Inhom,S,tm;
  if not lhs(sol)::function(name) then error "This is not an output from dsolve" end if;
  t:=op(1,lhs(sol));
  cs:=indets(sol,'suffixed'(_C,posint));
  if cs={} then return sol end if;
  Hom,Inhom:=selectremove(has,rhs(sol),cs);
  if not type(Inhom,`+`) then return sol end if;
  S:={};
  for tm in Inhom do
    if solve('identity'(Hom=tm,t))<>NULL then S:=S union {tm} end if
  end do;
  Inhom:=subs(S=~0,Inhom);
  lhs(sol)=Hom+Inhom
end proc:

You seem to have answered your own question.
Try this:

Student:-NumericalAnalysis:-Quadrature(sin(x), x=0..Pi, method=newtoncotes[open,4], output=information);

There is a discrepancy, however, between what is written in the help page
?Student:-NumericalAnalysis:-Quadrature
under Options and what is actually the case.
We read:
"adaptive = true or false
  Whether to apply the adaptive version of the specified method. The following methods are available when adaptive is set to true:
Newton-Cotes Rule, open or closed, with a user-specified order"

but if you try

Student:-NumericalAnalysis:-Quadrature(sin(x), x=0..Pi, adaptive=true,method=newtoncotes[open,4], output=information);

then you will get the error message

Error, (in Student:-NumericalAnalysis:-Quadrature) the newtoncotes[open, 4] method is not available when the adaptive option is specified

##
Note: I have submitted an SCR.

I changed pi into Pi and made an assignment to K in your worksheet.
I chose as an example to make all parameters in the system equal to 1 and all initial values equal to zero.
Change that as appropriate.
MaplePrimes18-03-07odesys.mw
To get the plots of all 4 separately you can use this code:

PLTS:=[seq(plots:-odeplot(res,[t,F(t)],0..7,numpoints=10000),F=[V1,V2,V3,V4])]:
plots:-display(Array(1..4,1..1,PLTS));

The plot at the end showing V1 and V2 only (and together):

You don't have a differential equation, but just two expressions u and m (where m happens to be the derivative of u w.r.t. t).
Just use a parametric plot.
 

restart;
u:=2*cos(3*Pi/2+4*t)/sqrt(64+(64-64*v^2)/v^2*exp(-eps*t));
v:=2: eps:=0.1:
simplify(u);
m:=diff(u,t);
plot([u,m,t=-2..10],labels=["u","u'"],numpoints=1000,color=gray);

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