restart;
S := Matrix( [ [a_1, a_2, a_3, a_4], [b_1, b_2, b_3, b_4], [c_1, c_2, c_3, c_4], [d_1, d_2, d_3, d_4] ] );
eq1:=a_1*d_1 = b_1*c_1, a_2*d_2 = b_2*c_2, a_3*d_3 = b_3*c_3, a_4*d_4 = b_4*c_4;
R := Matrix( [ [s_1, t_1, r_1, l_1], [s_2, t_2, r_2, l_2], [s_3, t_3, r_3, l_3], [s_4, t_4, r_4, l_4] ] );
eq2:=s_1*l_1 = t_1*r_1,s_2*l_2 = t_2*r_2,s_3*l_3 = t_3*r_3,s_4*l_4 = t_4*r_4;
T := Matrix( [ [1, 0, 0, 1], [0, 1, 0, 0], [0, 1, 0, 0], [1, 0, 0, 1] ] );
S.R=~T;
convert(%,set);
res:=solve(% union {eq1,eq2});
nops([res]);
res[1]; # First of the ninety

DE := diff(x(t), t) = x(t);                   
DF := diff(y(t), t) = k*y(t); #You missed k

xres:=dsolve(DE);
yres:=subs(_C1=_C2,dsolve(DF));
eq:=abs(y(t))=C*abs(x(t))^k;
eval(eq,{xres,yres}) assuming real;
combine(expand(%)) assuming real;

So indeed there is a constant C s.t. eq is satisfied (if that was what you wanted to show?).

Since the body of a procedure (function) is not evaluated until it is invoked, the loop variable i in your phi functions is just 'i' untill you actually make a call to phi. At that time i evaluates to the global value it has, which will be n+1 in your case (unless i is changed in other ways than by the loop).
So use unapply:
for i from 1 to 5 do phi[i]:= unapply((i+2)/(i+1)*x^j+x^(j+1),x) end do;

Try

indets(ratio,function);

and you will see expressions like

B[1](m[10]^2+m[20]^2)

These are of type function in Maple, same as f(x) or f(2).
The expression should have been
B[1]*(m[10]^2+m[20]^2)

You may try to repair the damage the easy way:

newratio:=evalindets(ratio,function,x->op(0,x)*op(1,x)):
eval(newratio, [m[10] = 0, m[20] = 0, m[30] = 1]);

You can use eval instead. Compare:

u:=sin(x)+cos(y);
subs({x=Pi/4,y=Pi/2},u);
%;
# with
eval(u,{x=Pi/4,y=Pi/2});

# and see

?subs
where it says:

" Simple applications of replacing a symbol by a value in a formula should normally be done with the eval command."

You need a parenthesis:
with(DEtools):

Since g(k)=Int(x(t)^2,t=0..k) we have that g'(k) = x(k)^2. Include that ode with g(0)=0:

restart;
eqn := x(t)^2+cos(x(t))+diff(x(t), t)+diff(x(t), t, t) = 100;
eqn2:=diff(g(t),t)=x(t)^2;
ics := x(0) = 2, (D(x))(0) = 3;
sol := dsolve({eqn, eqn2, ics,g(0)=0}, numeric, maxfun = 500000, output = listprocedure);
plots:-odeplot(sol,[t,g(t)],0..10);
#The procedures for x and g if needed:
X,G:=op(subs(sol,[x(t),g(t)]));

Using Maple 16 I get these 4 solutions with the same code:
1, 1/2, 1/4*(sqrt(5)+1)^2,  0;

Since Ts and Ths are lists of matrices and Q is a matrix, you could do
Qb:=[seq(Ts[i]^(-1).Q.Ths[i],i=1..nops(Ts))];
to get a list of the matrices you need.

Another way:

Qb:=Ts^~(-1).~[Q$nops(Ts)].~Ths;

I find the first one easier to understand.

The way Maple works this is not possible as you entered the commands.
No problem if you defined c before a and b or if later you defined c:='a+b'.

restart;
a:=5*x^2-3*y^3:
b:=3*x^2+10*y^3:
c:='a+b';
c;
eval(c,1);
restart;
c:=a+b;
a:=5*x^2-3*y^3:
b:=3*x^2+10*y^3:
c;
eval(c,1);

A simple attempt:

f:=(x,y)->x^2*sin(y);
g:=`if`(x+y<0 and x>0,undefined,f(x,y));
plot3d(f(x,y),x=-1..1,y=-Pi..Pi,axes=boxed);
plot3d(g,x=-1..1,y=-Pi..Pi,axes=boxed);

As you notice the one edge is jagged, the other straight. So this does not produce a perfect picture. But you could use a finer grid.

Edited: Notes in bold inserted.

I tried the following. It seems sound (No!) for even values of n; it is, however, way too slow for somewhat large n.
More importantly, it should only work for n even. (It should not work for any n).
But a first version surprisingly also "works" for odd n, i.e. gives the correct answer, but how?

restart;
with(IntegrationTools):
A:=n->Int(mul(sin(x/2^k)*2^k/x, k = 0 .. n), x = 0 .. infinity);

N:=10:
Change(A(N),x=2^N*t):
B:=combine(%):
f1:=GetIntegrand(B):
expand(f1,sin,cos): #Early version (may be commented out)
map(int,%,t=0..infinity); #Early version (may be commented out)
## In fact also in this case you had a sum of divergent integrals!
q:=op(1,f1);
f:=expand(f1/q,sin,cos):
Q:=Int(sin(n*t)/t^k,t=0..infinity);
Change(Q,t=u/n,[u]) assuming n>0;
Qres:=value(%) assuming k::posint,k>1;
##The relevant assumption should have been k>0,k <2
##The result also comes without assumptions.
S1:=map(Int,f,t=0..infinity): #A sum of divergent integrals!
S:=map(expand,S1,sin,cos):
eval(S,{seq( eval(Q=Qres,k=N+1),n=1..2^(N+1))})*q;
###################################################
N:=3:
Change(A(N),x=2^N*t):
B:=combine(%):
f1:=GetIntegrand(B):
expand(f1,sin,cos);
map(int,%,t=0..infinity);
#The integrals in the sum are all divergent!!
#To explore this further:
convert(expand(f1,sin,cos),list);
map(int,%,t=0..infinity);
convert(%,`+`);
#But directly:
int(cos(3*t)/t^4,t=0..infinity);

I have given a sound (hopefully) version in a comment below.

In your uploaded file For_prime.mw you provided an image. It would have spared me some time if you had typed in the integral.
However, I typed it.
I have uploaded a worksheet:

MaplePrimes13-02-16I.mw

You can do several simulations by doing a loop:

#Here I do 100.
for i from 2 to 101 do
     res := dsolve({ode2, ic[1]}, numeric, known = r(t), method = classical[foreuler], stepsize = 0.1e-1);
     p[i]:= odeplot(res, [[t, U(t)]], t = 365*k .. 365*k+365, colour = red);
end do:
display(seq(p[i],i=1..101));

And for fun:

panim:=display(seq(p[i],i=2..101),insequence=true):
display(p[1],panim);
#Will draw 3 curves: One corresponding to the mean value of r the other two deviating there from by 2*s to each side.
m:=-0.2e-1;
s:= 0.4e-1;
#Since I'm going to do 3 computations quite alike I use the parameters option:
odepar := diff(U(t), t) = -(A+q+B*U(t))*U(t);
respar:=dsolve({odepar, ic[1]}, range = 365*k .. 365*k+365, numeric,parameters=[q]);
respar(parameters=[m]);
q1:=odeplot(respar, [[t, U(t)]], t = 365*k .. 365*k+365, colour = green):
respar(parameters=[m+2*s*m]);
q2:=odeplot(respar, [[t, U(t)]], t = 365*k .. 365*k+365, colour = black):
respar(parameters=[m-2*s*m]);
q3:=odeplot(respar, [[t, U(t)]], t = 365*k .. 365*k+365, colour = black):
display(q1,q2,q3,panim,view=[90..100,170..230]);
display(q1,q2,q3,panim);
display(q1,q2,q3,seq(p[i],i=1..101));

#Whether this black band has anything to do with confidence = 0.95 I have no idea.
Is the question: "Is the chance 95% that the curves stay within the band all the time?"
or is the question: "Is the chance 95% that each curve stays inside the band 95% ot the time? "

Don't assign to xi(t), eta(t) and their derivatives. Instead do
eq1:= xi(t) = 1/2 - mu + ... ;
eq2:= eta(t) =  sqrt( ..... ;
and likewise with their derivatives.

Then do

eval(A,{eq1,eq2,eq3,eq4});