In order to get to the essential problem I simplied your equation some. The homogeneous boundary and initial conditions are the same.
My suggested solution is pde4 with res4 and the final plot.

restart;
pde1:=diff(w(x,t),t,t)+diff(w(x,t),t)+diff(w(x, t), x$4)=sin(t); #The simplified pde
ibc:={w(0, t) = 0, w(1, t) = 0, D[1, 1](w)(0, t) = 0, D[1, 1](w)(1, t) = 0,w(x, 0) = 0, D[2](w)(x, 0) = 0};
res1:=pdsolve(pde1,ibc,numeric);
res1:-plot(w(x,t),x=1/2,t=0..3*Pi);
pde2:={diff(w(x,t),t,t)+v(x,t)+diff(w(x, t), x$4)=sin(t),diff(w(x,t),t)=v(x,t)};
res2:=pdsolve(pde2,ibc,numeric);
res2:-plot([w(x,t),[v(x,t),color=blue]],x=1/2,t=0..3*Pi);
#pde3:={diff(w(x,t),t)+a(x,t)+diff(w(x, t), x$4)=sin(t),diff(w(x,t),t,t)=a(x,t)}; #As in your case this won't work
#res3:=pdsolve(pde3,ibc,numeric); #no good
#res3:-plot([w(x,t),[a(x,t),color=blue]],x=1/2,t=0..3*Pi); #no good
pde4:={diff(w(x,t),t,t)+v(x,t)+b(x,t)=sin(t),diff(w(x,t),t)=v(x,t),diff(w(x, t), x$4)=b(x,t)};
res4:=pdsolve(pde4,ibc,numeric);
#The factor 25 is just for scaling:
res4:-plot([25*w(x,t),[25*v(x,t),color=blue],[b(x,t),color=green]],x=1/2,t=0..3*Pi);
#According to pde4 the following is a plot of diff(w(x,t),t,t):
res4:-plot(sin(t)-v(x,t)-b(x,t),x=1/2,t=0..3*Pi);

There is a good reason for the difference. For your system the solutions for u, x, and z grow exponentially. Thus a tiny difference in initial conditions at t=0 from the correct one is likely to cause quite a difference at t = 4.
In your case the system can be solved exactly, i.e. we can find analytical expressions for the solutions. Thus we can see what goes on:

sys := {diff(u(t), t) = 4*u(t), diff(x(t), t) = 6*x(t)+2, diff(z(t), t) = x(t)+u(t)};
backw:={x(4)=10,u(4)=1,z(4)=2}; #Initial conditions at t = 0.
res1:=dsolve(sys union backw,numeric);
res1(0); #These will be used for solving in the forward direction:
forw:=subs(%,{u(0)=u(t),x(0)=x(t),z(0)=z(t)});
res2:=dsolve(sys union forw,numeric);
res2(4);
#output roughly [t = 4., u(t) = .86, x(t) = .16, z(t) = .33], which is way off as you say
#Now solve exactly:
res:=dsolve(sys union backw);#You notice the exponential growth
evalf(eval(res,t=0)); #Compare with:
forw; #Not at all bad, but different, and different enough to make a difference, when going forward
U,X,Z:=op(subs(res,[u(t),x(t),z(t)]));
plots:-odeplot(res1,[t,x(t)-X],0..4);
plots:-odeplot(res2,[t,x(t)-X],0..4,color=blue);

f:=x^5+a*x^4+b*x^3+c*x^2+d*x+e;
simplify(f,{x^3=1});
#See:
?simplify,siderels

This surprised me, maybe mainly because I rarely use Diff, but almost always diff.
I would have expected an answer like the one you get by using the other inert form for diff, viz. %diff:

restart;
test1 := (3*f-f^2)*(Diff(g, x))+(-f+3*f^2-2)*(Diff(g, y))+(2-f^2)*(Diff(g, z));
solve(test1=0,f); #A weird answer! (A bug maybe in solve?)
#However, %diff performs fine:
test2:=subs(Diff=%diff,test1);
solve(test2=0,f);
#################
#Or try substituting any unassigned name for Diff. Here I use G:
test3:=subs(Diff=G,test1);
solve(test3=0,f);
#You can get back to what you really want by substituting again:
subs(G=Diff,[%]);

## I filed an SCR
##And where is g, you ask. I think solve factors out g from Diff and gets
Diff(g,x) = g*Diff(1, x)
Then the answer is understandable, but hardly expected.

You need to use ` ` instead of ' ' as shown here:

restart;
`u/t`:=7; # `u/t` is a name and can be assigned to
'u/t':=67; # 'u/t' is not a name but is u/t surrounded by unevaluation quotes

Introduce a new unknown function u(y,zeta) = Theta(y,zeta) + mu*y^2/2 + mu*zeta, then u still satisfies the heat equation, but its derivatives w.r.t. y are zero at y=0 and y=1.
Then solve in textbook manner by separating variables.
pde:=diff(Theta(y,zeta),zeta)=diff(Theta(y,zeta),y,y);
ibc:={D[1](Theta)(0,zeta)=0,D[1](Theta)(1,zeta)=-mu, Theta(y,0)=1};
eval(pde,Theta(y,zeta)=u(y,zeta)-mu*y^2/2-mu*zeta);
pde2:=%+(mu=mu);
pdsolve(pde2);
ibc2:=eval(ibc,Theta=unapply(u(y,zeta)-mu*y^2/2-mu*zeta,y,zeta));

There is no way of getting an analytic expression using the usual mathematical functions. It can easily be solved for concrete values of C3, however. Obviously s = 0 is always a solution. Assuming that your C3 is real it is clear that there are no other solutions if C3 >= 1 or C3 <= 0.
For concrete values (here C3 = 1/2) you can do:
restart;
equ:= tanh(s)=C3*s;
C3:=1/2:
solve(equ,s);
allvalues(%);
evalf(%);
#Alternatively, use
fsolve(equ,s=2);
#or
RootFinding:-Analytic(equ,s=-3-I..3+I);

If you replace BCs by
BCs := {u(0,t)=sin(t), u(10,t)=0,T(0,t)=tanh(100*t^2), T(10,t)=0,u(x,0)=0,T(x,0)=0};
then T is continuous and even differentiable along the boundary. 
That makes quite a difference since you are interested in the behavior for x = 0.
It is not too surprising that discontinuous boundary conditions give rise to numerical problems.

Although not being familiar with the Physics package I tried this:

restart;
with(Physics):
with(Vectors):
_i;
type(_i,PhysicsVectors);
showstat(`type/PhysicsVectors`);

You are probably not thinking of true and false.
However, to me the following behavior was surprising.
There is automatic simplification in both cases, and it looks as if p1 and p2 simplify to the same procedure, but they don't.
restart;
p1:=proc(true) not false end proc;
p1(8); #returns true
showstat(p1);
p2:=proc(true) if not false then true else 67 end if end proc;
p2(8); #returns 8
showstat(p2);
#Similar constructions with parameter false behave similarly:
p3:=proc(false) not true end proc;
p4:=proc(false) if true then false else 67 end if end proc;

#More confusion:
restart;
p1:=proc(true) not false end proc;
p1(8); #Returns true
if % then 89 end if; #Shows that true is the boolean true
#Now try debugging
debug(p1);
p1(8); #Now it returns 8
undebug(p1);
p1(8); #It still returns 8
p1(78); #Returns 78
#Apparently debugging changes the address. If you try
addressof(eval(p1));
#before and after debugging you find different addresses, but that goes for p2 as well.
#More
#Making a copy of p1 before debugging:
q1:=eval(p1);
addressof(eval(q1)); #Same as for p1
#Debugging p1 doesn't alter the address of eval(q1). And
evalb(eval(p1)=eval(q1));
#returns false.
#However, if you do the similar things with p2 (with copy q2) then again after debugging p2, eval(p2) and eval(q2) have different addresses, but this time
evalb(eval(p1)=eval(q1));
#returns true, which maybe doesn't contradict what the help page for evalb says:
"The evalb command uses address tests to determine equality".
because this doesn't say that evalb doesn't use other tests too.
#Exhibiting the explicit difference
#Maybe the following throws some light on this rather confusing matter:
restart;
p1:=proc(true) not false end proc;
dismantle(eval(p1));
a:=addressof(eval(p1));
d:=disassemble(a);
debug(p1): undebug(p1):
a1:=addressof(eval(p1));
d1:=disassemble(a1);
evalb~([d]=~[d1]);
pointto(d[6]);
pointto(d1[6]);
dismantle(eval(p1));
#ToInert(eval(p1)) shows the same difference.

If A,B, and C are real at least one of the roots is certainly real. However, it is not necessarily the one appearing with no explicit 'I' in the result of solve:
p:=x^3+A*x^2+B*x+C;
res:=solve(p=0,x);
with(plots):
complexplot(eval([res],{B=1,C=1}),A=-2..2,style=point);
animate(complexplot,[eval([res],{B=1,C=1}),A=-2..A1,style=point],A1=-2.1..2);
animate(complexplot,[eval(res[1],{B=1,C=1}),A=-2..A1,style=point],A1=-2.1..2);


#The condition for exactly one real root is:
discrim(p,x)<0;
#which returns
-4*A^3*C+A^2*B^2+18*A*B*C-4*B^3-27*C^2 < 0
#############
To deduce that result in Maple you can use the usual method of changing the variable so that the coefficient of the second order term disappears:

p1:=subs(x=y-v,p);
collect(p1,y);
p2:=subs(v=A/3,%);
#Thus we consider:
p3:=y->y^3+a*y+b;
diff(p3(y),y);
#If a >= 0 then certainly only one real root since p3 is increasing.
#If  a < 0 then we demand that p3 be positive at the local minimum:
p3(sqrt(-a/3)) > 0;
simplify(%) assuming a<0;
#This may be written as
27*b^2+4*a^3>0;
#and has the advantage that it covers the case a >=0.
#So for our polynomial p the requirement for exactly one real root and two imaginary is:
subs({a=coeff(p2,y,1),b=coeff(p2,y,0)},%);
expand(%);
#Compare with
discrim(p,x)<0;

It helps to split the sum:

A:=Sum(I^n*BesselJ(n,2*sqrt(34))*exp(-I*n*arctan(5/3))*exp(I*n*Pi/4),n=1..infinity);
evalf(A); #No answer
#Splitting at n=N and exhibiting the two contributions:
N:=20:
evalf[15]([subsop(2=(n=1..N),A),subsop(2=(n=N+1..infinity),A)]);
#The sum
convert(%,`+`);


###############
For plotting purposes you can safely truncate the infinite series:
restart;
u:=BesselJ(n,r*sqrt(34))*exp(I*n*(Pi-arctan(5/3)+theta));
N:=40: #You may want to investigate how large it has to be (see below)
a:=unapply(Sum(u,n=1..N),r,theta);
applyop(evalc@abs,1,a(r,theta)); #Not needed but abs(a(r,theta)) is dominated by that.
plot3d(abs(a(r,theta)),r=0..1,theta=0..Pi/2);
###############
Using Abramowitz and Stegun, (9.1.20) we find a bound for BesselJ(n,R):
#Defining G to be
G:=(n,R)->R^n/2^n/GAMMA(n+1/2)*sqrt(Pi); #Can be improved to 2^(-n)*R^n/n!
#it follows from 9.1.20 that
G(n,R)-abs(BesselJ(n,R))>=0;
#assuming that n >=0 and R>=0.
#The tail end of your sum is dominated by
B:=N->Sum(G(n,r*sqrt(34)),n=N+1..infinity);
evalf(eval(B(20),r=1));
plot(B(20),r=0..1);
#The maximal value of r is very important. I have chosen 1. (you had infinity!).

Assuming that your answer is:
answer:=3*(1+q),
           (-3*m+sqrt(m*(256*m^3+160*m^2-31*m-16)))/(4*m*(m+1)),
           (-3*m-sqrt(m*(256*m^3+160*m^2-31*m-16)))/(4*m*(m+1));
then it is correct and agrees with Maple for q = 0, but not otherwise.
Example:
evalf(Eigenvalues(eval(M,{q=1,m=1})));
#Numerical computation of eigenvalues:
Eigenvalues(evalf(eval(M,{q=1,m=1}))); # Agrees with the above
evalf(eval(<answer>,{q=1,m=1})); #Doesn't agree with above
It is also somewhat suspicious that the first eigenvalue should depend on q only and the two other on m only.
Finally, you may try:
pol:=CharacteristicPolynomial(eval(M,m=1),lambda);
factor(eval(pol,lambda=answer[1]));

You could try solving numerically, but that would require concrete values for the parameters and you would have to specify initial conditions in addition to the boundary conditions.
You can only solve forward in time (think of the heat equation)!

restart;
sys2 := -diff(u(y, t), y, y) + S*diff(u(y, t), y) + diff(u(y, t), t) + M*u(y,t) + u(y,t)/k-theta(y,t) = 0,
 -diff(theta(y, t), y, y)/Pr + diff(theta(y, t), t) + S*diff(theta(y, t), y) = 0;
bcs1:=u(0,t)=0, theta(0,t)=1, u(1,t)=0, theta(1,t)=0; #As you have them for t>0.
inits:=u(y,0)=0,theta(y,0)=1; #Example
params:={S=1,M=1,k=1,Pr=2}; #Example
res1:=pdsolve(eval({sys2},params),{bcs1,inits},numeric,time=t,timestep=.001);
res1:-animate([[u,color=blue],[theta,color=red]],t=.4,0..1,frames=50,title="t=%f" );
#Now the computer doesn't explode if you try solving backwards, but almost:
bcs2:=u(0,t)=0, theta(0,t)=1, u(1,t)=1, theta(1,t)=0; # t<0
#Using the same initial conditions given at t = 0 as above:
res2:=pdsolve(eval({sys2},params),{bcs2,inits},numeric,time=t,timestep=.001);
res2:-animate([[u,color=blue],[theta,color=red]],t=-.4,0..1,frames=50,title="t=%f" );

You could try numerical solution (e.g. with the parameters option for the initial values at t0):

X:=Vector([seq(x[i](t),i=1..(L+1)^2)]);
diff~(X,t)=mat.X; #Your system
eqs:=convert(Equate(lhs(%),rhs(%)),set);
t0:=1: #Initial t
inits:={seq(x[i](t0)=x0[i],i=1..(L+1)^2)};
res:=dsolve(eqs union inits,numeric,parameters=[seq(x0[i],i=1..(L+1)^2)],range=1..15);
res(parameters=[1,2,3,.1,.2,.3,4,5,6]);
plots:-odeplot(res,[Re(x[1](t)),Im(x[1](t))],1..15,refine=1);