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

The first problem is that you begin by assigning to mu after which you define the equations. Thereby the equations won't change if you change mu later.
The next problem is that even if you postpone assigning to mu till you use dsove/numeric then that sol will be the solution corresponding to the assigned value of mu.

The solution is to use the parameters option in dsolve/numeric, see
?dsolve,numeric,interactive

So with the necessary changes you got:

restart;
with(plots):
with(DEtools):
k := 100; L := .25; M := .5; g := 4;
F[N] := k*L*(1-L/sqrt(x(t)^2+M^2));
eq1 := diff(x(t), t) = v(t);
eq2 := diff(v(t), t) = k*x(t)*(1-L/sqrt(x(t)^2+M^2))+mu*abs(F[N]-g)*signum(x(t));
sol := dsolve({eq1, eq2, v(0) = 1, x(0) = .1098},numeric,parameters=[mu]);
sol(parameters=[0.2e-1]); #Here mu will be replaced by that value
newplots := [odeplot(sol, [t, x(t)], t = 0 .. 2, colour = red), odeplot(sol, [x(t), v(t)], t = 0 .. 2, colour = blue)];
display(array(newplots));
#Now your loop:
for i to 17 do
    mu := (i-9)*(1/2); #now assignment is harmless
    sol(parameters=[mu]);
    sss := cat("The value of &mu; is ", convert(mu, string), ".");
    P[i] := odeplot(sol, [t, x(t)], t = 0 .. 2, colour = red, title = sss)
end do:

display(seq(P[k], k = 1 .. 17), insequence = true);

Using equation 3.4.82 in the reference you gave and also
diff(phi,t,t) = phidot*diff(phidot,phi) we are back to the second order ode:

ode:=diff(phi(t),t,t)=-(sqrt(3/2)/M^2*sqrt(diff(phi(t),t)^2+m^2*phi(t)^2)*diff(phi(t),t)+m^2*phi(t));
#Following Adri van der Meer we turn the equation into a first order system:
sys:=diff(phi(t),t)=phi1(t),subs(diff(phi(t),t)=phi1(t),ode);
#I pick my own parameters, m = M = 1:
with(DEtools):
DEplot(eval([sys],{M=1,m=1}),[phi(t),phi1(t)],t=0..30,[[phi(0)=1,phi1(0)=0],[phi(0)=-1,phi1(0)=0],[phi(0)=2,phi1(0)=-1],[phi(0)=-2,phi1(0)=1]],color=grey,linecolor=blue,stepsize=.1,thickness=1);