@mehran1520 In equations EQ4 and EQ5 there appears the parameter T0. It needs to have a numerical value. In my example above I used T=1, since I had no way of knowing its value. But I suppose you do?

You cannot specified the derivatives for EQ4 and EQ5. In those equations the derivatives of highest order appearing are of order 1 and are of q4 and q5 respectively. So you must limit yourself to
inits:={q1(0) = 0.1e-2, q2(0) = 0, q3(0) = 0, q4(0) = 0, q5(0) = 0, D(q1)(0) = 0, D(q2)(0) = 0, D(q3)(0) = 0};

@bastorer Maybe using copy~ like this should suffice:

restart;
Array1 := Array([Array([1, 2]), Array([Array([8,9]), 4])]); #Example nested twice
Array2:=((copy~)~)~(Array1); #Overdone once, but doesn't hurt!
Array1[1][1]:= 11:
Array1;
Array2;
Array1[2][1][2]:=Pi;
Array1;
Array2;

It is interesting that the result of evalf(dsolve({sys})) only contains 3 arbitrary constants. There ought to be 4.

The system can be rewritten as a linear first order system of the form X'(x) = M.X(x) + C.
Its solution can be expressed in terms of eigenvectors and eigenvalues of M which are easily found and are manageable since the coefficients of M are floats.
Doing that I found a solution which agrees reasonably well with the solution obtained by Carl.

restart;
A := 54.15836673*(diff(X2s(x), x, x)) = -365.4395362*(diff(X1s(x), x))+208.2315661*X2s(x);
B := 641.1196154*(diff(X1s(x), x, x)) = 365.4395362*(diff(X2s(x), x))-2.575699975*X1s(x)-7.882173342;
bc := X1s(0) = 0, X1s(15) = 0, X2s(0) = 0, X2s(15) = 0;
sys := A, B:
dsolve({sys}):
indets(%,name);
cvs:=DEtools[convertsys]([sys],{},[X1s(x),X2s(x)],x,Y,Yp);
rhs~(cvs[1]);
M,C:=LinearAlgebra:-GenerateMatrix(rhs~(cvs[1]),[seq(Y[i],i=1..4)]);
C:=-C; #Notice
Lambda,V:=LinearAlgebra:-Eigenvectors(M):
Lambda:=simplify(fnormal~(Lambda));
V:=simplify(fnormal~(V));
A particular solution to X' = M.X+ C is
Xp:=simplify(-M^(-1).C);
The general solution is
X(x)= V.<seq(exp(Lambda[i]*x)*c[i],i=1..4)>+Xp;
gensol:=rhs(%);
eq1:=eval(gensol,x=0);
eq2:=eval(gensol,x=15);
solve({eq1[1]=0,eq1[3]=0,eq2[1]=0,eq2[3]=0},{seq(c[i],i=1..4)});
sol:=eval(gensol,%);
So the solution is
X1s(x)=sol[1];
X2s(x)=sol[3];
odetest({X1s(x)=sol[1],X2s(x)=sol[3]},{sys});
fnormal(%,8);

I tried a reduced system choosing the first n such that seq(eq[i],i=1..n) is a nonlinear system when seq(a[i]=0,i=n+1..9) is used.

n needed to be 4.
Even here you have a singularity. A singularity is not surprising since there are several quadratic terms and even one cubic.

solve({seq(eq[i]=0,i=1..4)},{seq(diff(a[i](t),t,t),i=1..4)});
sys4:=eval(%,{seq(a[i]=0,i=5..9)});
res4:=dsolve(sys4 union op~({seq([a[i](0)=0,D(a[i])(0)=0],i=1..4)}),numeric);
odeplot(res4,[t,a[1](t)],0..0.001);

Following your description in the worksheets, you intend to solve the system

M.V22+K(t).V+P(t) = EQ

where 
EQ:=<seq(eq[i],i=1..9)>:
V22:=<seq(diff(a[i](t),t,t),i=1..9)>:
V:=<seq(a[i](t),i=1..9)>:
#and M, K(t), P(t) as described in the worksheets.
Now you say that you define K(t) such that
M.V22+K(t).V+P(t) = EQ holds.
and in fact you did:
simplify(M.V22+K(t).V+P(t) - EQ);
returns the zero vector.
But then what you are trying to solve according to your desription is equivalent to solving 0=0.
Thus I don't understand what is going on.

If you meant that K should be defined such that
M.V22+K(t).V+P(t) = 0 is equivalent to EQ=0,

then what you are trying to solve is the system EQ=0, i.e. in Maple
a:='a':
res:=dsolve(op~({seq([eq[i]=0,a[i](0)=0,D(a[i])(0)=0],i=1..9)}),numeric);
odeplot(res,[t,a[6](t)],0..0.4);
#The Maple solution becomes singular at roughly t=0.001.
I don't know anything about the Newmark method you mention.

@emmantop I don't get any error at all. I'm also using Maple 17.

Please try:

MaplePrimes14-03-17o.mw

@J4James Do you mean something like this:

restart:
Eq1:=diff(psi(y),y$4)-diff(psi(y),y$2)=0;
res1:=dsolve(Eq1);
bcs:=psi(h1)=F,D(psi)(h1)=-1,psi(h2)=-F,D(psi)(h2)=-1;
res2:=simplify(dsolve({Eq1,bcs},psi(y)));
match(rhs(res2)=rhs(res1),y,s);
s;

@Carl Love This is a somewhat weird system. I wonder how it originated.
You cannot isolate the highest derivatives diff(u(t),t$3) and diff(w(t),t$4).
Immediately from looking at it I thought well, w and u appear differentiated to highest order 4 and 3 respectively, so you need 3+4 = 7 conditions not 6.

This extremely simple example may illustrate the point:

ode1:=diff(x(t),t)=y(t);
ode2:=diff(x(t),t,t)=diff(y(t),t);
solve({ode1,ode2},{diff(x(t),t,t),diff(y(t),t)}); #No solution
dsolve({ode1,ode2}); #Quite a few
dsolve({ode1,ode2,x(0)=1,y(0)=1}); # None it says (silently
#However
dsolve({ode1,x(0)=1});
eval(ode2,%);

@Carl Love You are right. Originally (a long time ago) the dependent variable had to be given and in the form of (example) x(t), thereby revealing both x and t. Now it is optional, but necessary in the present case but also in a case like

ode:=diff(x(t),t)=diff(h(t),t)+x(t);
dsolve({ode,x(0)=0});
dsolve({ode,x(0)=0},x(t));

I would recommend using unapply as in

g:=unapply(L(f,f)(x),x);

This has the advantage that the integration takes place at this step, whereas using

g:= x->L(f,f)(x);

no integration takes place. It takes place when a call to g is made as in g(x).

If phi(1) is negative (for some values of the parameters) and you think that your problem demands that it be positive, then either
1) your code or model is deficient
or
2) there are multiple solutions to your boundary value problem, which certainly cannot be excluded.

@samiyare I'm assuming that this is the same comment as the one I just replied to.

@samiyare My comment about your replacement for infinity, which you called binfinitive, was just meant to point out that in the context its role is just like any other parameter and may cause convergence problems as the other parameters might. I didn't mean to suggest that changing variables as I did would help anything at all as far as convergence is concerned.

However, it does open up the possibility of using continuation in inf from value like 6 to N (e.g. 15) with continuation=6*(1-c)+c*N, where we replace inf with 6*(1-c)+c*N.

Something like this:

restart;
MUR:=(1-phi)^(5/2):
RhoUR:=(1-phi+phi*rho[p]/rho[f]):
RhoCPR:=(1-phi+phi*rhocp[p]/rhocp[f]):
BetaUR:=(phi*rho[p]*beta[p]+(1-phi)*rho[f]*beta[f])/(RhoUR*rho[f])/beta[f]:
dqu3:=diff(h(x),x$1)-RhoUR*BetaUR*T(x);
dqu2:=5*diff(T(x),x$2)+k[f]/k[nf]*Pr*RhoCPR*f(x)*diff(T(x),x$1);
dqu1:=5/(MUR)*diff(f(x),x$3)
+ 2*(diff(h(x),x$1)*x-h(x))
+RhoUR*(3*f(x)*diff(f(x),x$2)-diff(f(x),x$1)^2);
PDEtools:-dchange({f(x)=fy(y),T(x)=Ty(y),h(x)=hy(y),x=y*inf},[dqu1=0,dqu2=0,dqu3=0],[fy,Ty,hy,y]);
sysy:=solve(%,{diff(fy(y),y,y,y),diff(Ty(y),y,y),diff(hy(y),y)});
bcsy:={Ty(0)=1,Ty(1)=0,fy(0)=0,D(fy)(0)=inf*lambda,D(fy)(1)=0,hy(1)=0};
rho[f]:=998.2: cp[f]:=4182: k[f]:=0.597:   beta[f]:= 2.066/10000:
rho[p]:=3380: cp[p]:=773: k[p]:=36:   beta[p]:= 8.4/1000000:
k[nf]:=((k[p]+2*k[f])-2*phi*(k[f]-k[p]))/((k[p]+2*k[f])+phi*(k[f]-k[p])):
rhocp[nf]:=rho[p]*cp[p]*phi+rho[f]*cp[f]*(1-phi):
rhocp[p]:=rho[p]*cp[p]:
rhocp[f]:=rho[f]*cp[f]:
N:=15: #So here inf is really going to be 15 by continuation from 6 to 15.
sysN:=subs(inf=6*(1-c)+N*c,sysy);
bcsN:=subs(inf=6*(1-c)+N*c,bcsy);
phi:=0.0:
Pr:=7: lambda:=0.1:
res:=dsolve(sysN union bcsN,numeric,continuation=c,output=listprocedure);
[seq([y,1/N^i*diff(fy(y),[y$i])],i=0..2)];
plots:-odeplot(res,[seq([y,1/N^i*diff(fy(y),[y$i])],i=0..2)],0..1,tickmarks=[[seq(i=N*i,i=0..1,1/N)],default],labels=[x,"f, f',f''"]);
plots:-odeplot(res,[y,fy(y)-y],0..1,tickmarks=[[seq(i=N*i,i=0..1,1/N)],default],labels=[x,"f(x)-x"]);
#Alternatively:
F,F1,F2:=op(subs(res,[seq(diff(fy(y),[y$i]),i=0..2)]));
plot([F(x/N),F1(x/N)/N,F2(x/N)/10^2],x=0..15,labels=[x,"f, f', f''"]);
plot(F(x/N)-x/N,x=0..N);

@Carl Love But would you consider the non-working of Q1 a bug?

(I just submitted an SCR).
I noticed that if in Q1 the output procedure is equipped with $ then the elementwise version works:

Q11:=proc(Var::list,$) local n;
  n:=nops(Var);
  print(Var,n);
  proc(var::list,$)
       subs(Var=~[seq(1..n)],var); #Elementwise
       #subs(zip(`=`,Var,[seq(1..n)]),var);
  end proc;
end proc:
Q11([x,y]); # works

@Carl Love Thanks for pointing that out. Obviously, it should have been D(fy)(1)=0.
I have corrected it above.