Obviously you don't mean literally that the first six arguments to pdsolve should be the integers 1,2,3,..,6.
You must be referring to labels or whatever.
However, taking your equations as they are and adding these two (and forgetting eq3):
eq7:= s1(t)=8*(diff(x(t), t))*h(t)+16*y(t)*x(t)+8*u(t)*(diff(y(t), t));
eq8:=diff(s(t),t)=s1(t);
#I did the following:
res:=dsolve({eq2,eq4,eq5,eq6,eq7,eq8,h(0) = 0.1e-10, s(0) = .6, u(0) = .62, x(0) = 596561.43714, y(0) = .12096},numeric,output=listprocedure);
U,S,S1,X:=op(subs(res,[u(t),s(t),s1(t),x(t)]));
plots:-odeplot(res,[[t,x(t)],[t,5000*y(t)],[t,u(t)],[t,100*s(t)],[t,s1(t)]],0..0.457,view=-50..75);
#Notice that there seems to be a singularity at t = 0.45739102
#That was just solving the odes.
#Now for the pde, i.e. eq1:
pde:=eval(eq1,{s(t)=S(t),u(t)=U(t),diff(u(t),t)=X(t),diff(s(t),t)=S1(t)});
#You need a boundary condition for pdsolve,numeric to get going. I chose T(t,0)=0:
sol:=pdsolve(pde,{T(0,r)=0,T(t,0)=0},numeric,time=t,range=0..125,timestep=0.01,spacestep=.5);
sol:-plot(t=.4);
sol:-animate(t=0.45,view=-50..50);
#At least we don't get any errors!

The error is that the right hand sides of your odes are row vectors (of dimension one).
Try
whattype(rhs(sys12));
op(1,rhs(sys1));
You should be able to repair that by doing:
SYS:=map(x->(lhs(x)=rhs(x)[1]),{sys1,sys2,sys3,sys4,sys5,sys6,sys7,sys8,sys9,sys10,sys11,sys12});
FF2:=dsolve(SYS union {init},numeric);
#However, this doesn't seem to finish within my time tolerance.

In Maple solving pdes numerically can be done for problems with 2 variables only.
See
?pdsolve,numeric
where it describes the first argument to pdsolve PDEsys as
"single or set or list of time-dependent partial differential equations in two independent variables".

You may try symbolically:
pdsolve(pde);
eval(%,u=(x->x)); #The case u(x)=x
value(%);
In the output appears an arbitrary function _F1 of 2 variables.
The case u(x)=x-x^3 gives a result with the same structure, but the second argument to _F1 is rather complicated.
Now the job is to find _F1 so that your other requirements are satisfied.
The form of the answer when u(x)=x makes a change of variables to polar coordinates seem a good idea:
pde2:=PDEtools:-dchange({x=r*sin(theta),p=r*cos(theta),f(x,p,t)=g(r,theta,t)},eval(pde,u(x)=x),[r,theta,g]);
pdsolve(pde2);
From the initial condition _F1 is immediately found. The solution found is always positive, thus doesn't satisfy your zero boundary conditions.

You got mixed derivatives that are not equal:
sys:=[
Diff(H(x1,x2,lambda), x1) = a1*x1 + a3*lambda,
Diff(H(x1,x2,lambda), x2) = b2*x2 + b3*lambda,
Diff(H(x1,x2,lambda), lambda) = c1*x1 + c2*x2
];
#While the first two are in agreement:
diff(sys[1],x2);
diff(sys[2],x1);
# these are not (unless a3=c1):
diff(sys[1],lambda);
diff(sys[3],x1);
#There is the same problem with this pair:
diff(sys[2],lambda);
diff(sys[3],x2);
##With the necessary requirements satisfied you get a solution:
sys2:=eval(sys,{a3=c1,b3=c2});
pdsolve(sys2);

You can construct a function with a remember table by doing the following:

f:={[1,2],[2,a],[3,b]}; #You need lists []
assign(seq(F(p[1])=p[2],p=f));
F(1);
F(2);
F(3);
op(4,eval(F));

ode:=diff(x(t),t)=x(t);
res:=dsolve({ode,x(0)=1},numeric,output=Array([seq(.01*i,i=0..100)]));
A:=res[2,1]; #The matrix of t,x values
plot(A); #Testing
ExportMatrix("F:/testMatrix.txt",A);
#The file can be opened in any text editor

With Digits=10 (default) I get 3 solution sets for hal. With Digits=15 I get one set. With Digits=25 I get 3.
Try converting to rational:
Let your final equations be eqs. Then try
Eqs:=convert(eqs,rational);
Hal := solve(Eqs, {C0, C1, C2, k0, k1, k2});
nops([Hal]); #Answer 3 when Digits=10, (but 1 when Digits=15)
## Checking all 3 solutions:
seq(simplify((lhs-rhs)~(eval(Eqs,Hal[i]))),i=1..3);

You can include the equation for s(t) in the system. Say your system including initial conditions is 'sys'. Then do
sys2:={op(sys),s(t)=diff(u(t)*h(t),t)}; #I'm leaving the unknown r out
sol2:=dsolve(sys2,numeric);
plots:-odeplot(sol2,[t,s(t)],1.44..2.4);

Since there is no general formula for the roots of polynomials of order 5 and above Maple cannot (in general) do better than it did in your case.
You can try numerical solution instead using fsolve with the extra argument complex.

You can certainly find numerical solutions. You need initial (or boundary) conditions though.
Example:
res:=dsolve({DE1,DE2,Y(0)=1,b(0)=2},numeric);
plots:-odeplot(res,[[t,b(t)],[t,Y(t)]],-.3..0.5);
plots:-odeplot(res,[Y(t),b(t)],-.3..0.5);#Phase space plot
##Equilibrium point:
solve(rhs~({DE1,DE2}),{b(t),Y(t)});
allvalues(%);
evalf(%);
#Phase space plot of two orbits starting close to the equilibrium point:
DEtools[DEplot]({DE1,DE2},[Y(t),b(t)],t=-1..10,[[Y(0)=0.5,b(0)=1.4],[Y(0)=0.52,b(0)=1.4]],Y=0..2,b=0..2,stepsize=0.1);

To see what is going on it might be a good idea to use implicitplot.
First with your worksheet as it is:
implicitplot(h=0,a=0..1,P=0..4,gridrefine=2);
#Then after restart with the assignment to S commented out:
implicitplot(eval(h,S=0.5)=0,a=0..1,P=0..4,gridrefine=2);
animate(implicitplot,[h=0,a=0..1,P=0..4,gridrefine=2],S=0..1);

Root finding could be done with RootFinding:-NextZero like this, ehere the inputs to PC are the actual values of the unassigned global names S and a.

PC:=proc(s,aa) local k,r,N;
   N:=5;
   r[0]:=0;
   for k from 1 to N do
      r[k]:=RootFinding:-NextZero(unapply(eval(h,{S=s,a=aa}),P),r[k-1]);
      if r[k]=FAIL then N:=k-1; break end if
   end do;
   [seq(r[i],i=1..N)]
end proc;
PC(1,.2);
PC(0.5,.4);

I believe the following gives you the 3 lowest eigenfrequencies. Notice the number of wiggles for each one.

restart;
EOM:=diff(EI(x)*diff(phi(x),x,x),x,x)-rho*A(x)*lambda*phi(x)=0;
EI:=x->10^(-5)*x^3+10;
A:=x->x/80+3;
rho:=8000;
k_h:=9000000;
BC1:=convert(eval(diff(EI(x)*diff(phi(x),x,x),x)=k_h*phi(x),x=88),D);
BC2:=EI(88)*(D@@2)(phi)(88)=k_r*D(phi)(88);
BC3:=convert(eval(diff(EI(x)*diff(phi(x),x,x),x),x=0)=0,D);
BC4:=EI(0)*(D@@2)(phi)(0)=0;
dsys:={EOM,BC1,BC2,BC3,BC4};
res1:=dsolve(eval(dsys,{k_r=k_h/2}) union {phi(0)=100},numeric);
plots:-odeplot(res1,[x,phi(x)],0..88);
res1(0);
omega1:=sqrt(subs(res1(0),lambda));
res2:=dsolve(eval(dsys,{k_r=k_h/2}) union {phi(0)=100},numeric,approxsoln=[phi(x)=100,lambda=.0000001]);
plots:-odeplot(res2,[x,phi(x)],0..88);
omega2:=sqrt(subs(res2(0),lambda));
res3:=dsolve(eval(dsys,{k_r=k_h/2}) union {phi(0)=100},numeric,approxsoln=[phi(x)=100,lambda=.000000005]);
plots:-odeplot(res3,[x,phi(x)],0..88);
omega3:=sqrt(subs(res3(0),lambda));
omega1,omega3,omega2; #Increasing order

You probably didn't mean to have  eval(PSI,x=1.366)=1 twice, but meant the second to be
eval(PSI,x=-1.366)=1 (?).
Assuming that I would first try numerical solution:
a:=1.366:
res:=dsolve({VR22,VS22,phi(a)=1,D(phi)(a)=0,phi(-a)=1,D(phi)(-a)=0,psi(a)=1,psi(-a)=1},numeric);
plots:-odeplot(res,[[x,phi(x)],[x,psi(x)]],-a..a); p1:=%:
#After that continue as you did to get SOL with the correction I mentioned and then redefining PHI and PSI:
PHI,PSI:=op(subs(SOL,[phi(x),psi(x)]));
plot([PHI,PSI],x=-a..a);p2:=%:
plots:-display(p1,p2); #Indeed they fall right on top of each other.

sys:=solve({eq1,eq2,eq3,eq4},{diff(alpha(t), t, t),diff(theta(t), t, t),diff(x(t), t, t),diff(z(t), t, t)}):
subs(t=0,CI,denom~(rhs~(sys)));
%;
Initially you have alpha'' = undefined (and theta'', x'', z'') because of division by zero.
Start by trying other initial conditions, e.g. just changing theta(0)=0 to theta(0)=0.00001.

You have two problems.
1. Starting a worksheet with local declarations valid for that session was introduced (I believe) in Maple 17.
Certainly it is not possible in Maple 12.
2. The elementwise operation ~ was introduced in Maple 13. Instead you can use map.

You could take a look at overload:

?overload
restart;
MatrixOperations := module() option package;  export `+`;
    `+` := proc(a::Matrix, b::Matrix) option overload;
        :-`+`(map(x->x^2,a),map(x->x^2,b))
    end proc;
end module;
with(MatrixOperations);
A:=Matrix([[1,2,3],[3,4,5]]);
B:=Matrix([[11,22,33],[33,44,55]]);
A+B;