I don't get any such Z1 in Maple 17,16,15.

I haven't updated the Physics package. I get no such error in Maple 17.02.

@Stephan This seems to me to be a bug.
You can use isolate instead:

isolate(eq[1],sigma[F](xi));

To see that it is a bug you can try

restart;
p:=piecewise(x<1,67,0);
solve(p=a/b(0),a);
solve(p=a/b,a);
solve(p=a/exp(1),a);


I'm using Maple 17 as you are. I don't know if this is fixed in Maple 18.
The bug is NOT present in Maple 15 nor in Maple 16.

Just for the record I shall submit an SCR (bug report).

@J4James I think I finally got what you want:

restart;
ODE:=diff(u(y),y)=y*m*exp(Lambda*(y*m)^2);
bcs:=D(u)(0)=0,u(h)=-1;
dsolve({ODE,bcs});
uu:=rhs(%);
us:=simplify(F-int(uu,y=0..h)) assuming Lambda>0;
h:=1+phi*cos(2*Pi*x):
F:=A-1:
Lambda:=2*B*W^2:
A:=2:x:=1:B:=3:W:=2:
uu;
U1:=evalf(convert(us,dawson));
p:=proc(phi1) if not type(phi1,numeric) then return 'procname(_passed)' end if;
   fsolve(eval(U1,phi=phi1)=0,m)
end proc;
p(0);
p(0.2);
U:=phi1->evalf(eval(uu,{phi=phi1,m=p(phi1)}));
U(0);
plot([U(0),U(0.2)],y=-1..1,axes=box,color=[red,green]);
plots:-animate(plot,[U(phi),y=-1..1,axes=box],phi=0..0.5);

@J4James
restart;
ODE:=diff(u(y),y)=y*m*exp(Lambda*(y*m)^2);
bcs:=D(u)(0)=0,u(h)=-1;
dsolve({ODE,bcs});
uu:=rhs(%);
us:=simplify(F-int(uu,y=0..h)) assuming Lambda>0;
h:=1+phi*cos(2*Pi*x):
F:=A-1:
Lambda:=2*B*W^2:
A:=2:x:=1:B:=3:W:=2:
U:=evalf(convert(us,dawson));
p:=proc(phi1) fsolve(eval(U,phi=phi1)=0,m) end proc;
p(0);
p(0.2);
plot(p,0..0.2,adaptive=false,numpoints=5);

But I don't understand where your variable y is coming from in your "goal".

@mehdi jafari If you change the output in your procedure to

Array(1..Fx,1..Fy,1..Ft,(i,j,k)->u[i,j,k],datatype=float)

then you can use matrixplot as in

plots:-matrixplot(Matrix(M[..,..,5]),labels=[x,y,u],caption=typeset("Time fixed at t =%1",5));

And an animation in time:

S:=seq(plots:-matrixplot(Matrix(M[..,..,i]),labels=[x,y,u],caption=typeset("Time t =%1",i)),i=1..21):
plots:-display(S,insequence);

Notice that time is just represented by the index.

@J4James I don't think that there is any way of getting a solution for m in closed form.
You need values for all parameters to be able to find the (unique) solution there is for m.

@Stephan Yes you can easily remove them using evalindets. But how did they get there in the first place?
It might be better to correct the problem at that stage.
However, here is the removal by evalindets:

evalindets(FuncA,list,op);

@Stephan 

Funct:=piecewise(xi*L<=L/2,x^2,L/2>xi*L,x);
Funct2:=eval(Funct,L=1); #The extremely simple solution for your specific example
evalindets(Funct,relation,s->map(x->x/L,s)); #the more complicated
evalindets(Funct,`<=`,s->map(x->x/L,s)); #restricted to type `<=`

A slightly more complicated example shows the broader scope of using evalindets:

G:=piecewise(xi*L<=L/2,x^2+L*x,L/2>xi*L,x+sin(x*L));
evalindets(G,relation,s->map(x->x/L,s));

@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