You could try:

restart;
local `+`;
`+`:=proc(a,b) :-`+`(a^~2,b^~2) end proc;
1+2;
A:=Matrix([[1,2,3],[3,4,5]]);
B:=Matrix([[11,22,33],[33,44,55]]);
A+B;

Something like this skeleton:

Initialize := proc (p, theta, PiV, PiU, n, m) local J,A;
  J := CompJ(PiU, m);
  ####
  A:=[7,9,13];
  ###
  J,A
end proc;

Initialize(1,2,3,4,5,6);
#Assigning either like this
res:=Initialize(1,2,3,4,5,6);
res[1];
#or like that:
res1,res2:=Initialize(1,2,3,4,5,6);
res2;

It appears that there are several solutions.
Consider the following two rather different solutions obtained by using approxsoln.

restart;
ode:=
     diff(f(eta),eta$3) + f(eta)*diff(f(eta),eta$2) - (diff(f(eta),eta))^2 +
     lambda*(2*f(eta)*(diff(f(eta),eta))*(diff(f(eta),eta,eta)) -
            f(eta)^2*diff(f(eta),eta,eta,eta)) -
     k*(diff(f(eta),eta)-1) =
     0;
params:= [lambda= 0.1, k= 0.2, inf= 10]:
BCs:= f(0) = 0, D(f)(0) = 1, D(f)(inf) = 0:
asol:=2.5*sin(eta*Pi/2/10): #Picked to look like Carl's solution
res1:= dsolve(eval({ode, BCs}, params), numeric, approxsoln=[f(eta)=asol]);
plots:-odeplot(res1, [eta, f(eta)], 0..10);
#Now give the derivatives explicitly as 0 and 0:
res2:= dsolve(eval({ode, BCs}, params), numeric, approxsoln=[f(eta)=asol,diff(f(eta),eta)=0,diff(f(eta),eta,eta)=0]);
plots:-odeplot(res2, [eta, f(eta)], 0..10);
#My guess is that there are infinitely many solutions.

With maxfun=10^7 and compile=true the code runs smoothly.

restart;
with(plots):
mb:=765; mp:=587;Ib:=76.3*10^3;Ip:=7.3*10^3; l:=0.92; d:=10; F:=-1.2; omega:=0.43;g:=9.81;ly:=3;k:=0.02001014429;
A:=168913.8672;
s:=0.0666666666667;
n:=49.97465213;
eq1:=(mp+mb)*diff(x(t),t$2)+mp*(d*cos(theta(t))+l*cos(alpha(t)+theta(t)))*diff(theta(t),t$2)+mp*l*cos(alpha(t)+theta(t))*diff(alpha(t),t$2)+mp*(d*diff(theta(t),t)^2*sin(theta(t))+l*(diff(theta(t),t)+diff(alpha(t),t))^2*sin(alpha(t)+theta(t)))+A*2*(s*sinh(k*ly+k*ly)*sin(omega*t-k*x(t)))=0;

eq2:=(mp+mb)*diff(z(t),t$2)-mp*d*(sin(theta(t)+alpha(t))+sin(theta(t)))*diff(theta(t),t$2)-mp*l*sin(alpha(t)+theta(t))*diff(alpha(t),t$2)+mp*(d*diff(theta(t),t)^2*cos(theta(t))+l*(diff(theta(t),t)+diff(alpha(t),t))^2*cos(alpha(t)+theta(t)))+9.81*(mp+mb)+1000*g*z(t)*15.3*30+A*cosh(k*ly+k*z(t))*n*(cos(omega*t-k*15)-cos(omega*t+k*15))=0;
eq3:=mp*(d*cos(theta(t))+l*cos(alpha(t)+theta(t)))*diff(x(t),t$2)-mp*(l*sin(theta(t)+alpha(t))+d*sin(theta(t)))*diff(z(t),t$2)+(Ip+Ib+mp*(d^2+l^2)+2*mp*d*l*cos(alpha(t)))*diff(theta(t),t$2)+(Ip+mp*l^2+mp*d*l*cos(alpha(t)))*diff(alpha(t),t$2)-mp*sin(alpha(t))*(l*d*diff(alpha(t),t)^2-l*d*(diff(alpha(t),t)+diff(theta(t),t))^2)+mp*9.81*l*sin(alpha(t)+theta(t))+mp*9.81*d*sin(theta(t))=0;

eq4:=mp*l*cos(alpha(t)+theta(t))*diff(x(t),t$2)-mp*l*sin(alpha(t)+theta(t))*diff(z(t),t$2)+(Ip+mp*l^2+mp*d*l*cos(alpha(t)))*diff(theta(t),t$2)+(Ip+mp*l^2)*diff(alpha(t),t$2)-mp*9.81*l*sin(alpha(t)+theta(t))+l*d*mp*diff(theta(t),t$1)^2*sin(alpha(t))=0;
CI:= x(0)=0,z(0)=0,theta(0)=0,alpha(0)=0,D(x)(0)=0,D(alpha)(0)=0,D(z)(0)=0,D(theta)(0)=0;

solution:=dsolve([eq1,eq2,eq3,eq4,CI],numeric,maxfun=10^7,compile=true);

odeplot(solution,[[t,x(t)],[t,alpha(t)],[t,z(t)],[t,theta(t)]], t=0..1000, thickness=2);

odeplot(solution,[[t,x(t)]], t=0..1000, thickness=2);

odeplot(solution,[[t,z(t)]], t=0..1000, thickness=2);
odeplot(solution,[[t,alpha(t)]], t=0..1000, thickness=2);
odeplot(solution,[[t,theta(t)]], t=0..1000, thickness=2);

Apart from the constant solution n(t) = 1.036007094 I doubt that you are going to get an analytical solution.
You can try
dsolve(ode3);
but the result could hardly be called an analytical solution.
Finding constant solution(s):
eval(ode3,diff(n(t),t)=0);
solve(%,n(t));
#Solving numerically:
res:=dsolve({ode3,n(0)=1,D(n)(0)=0},numeric);
plots:-odeplot(res,[t,n(t)],-15..5);
#Plot in the phase plane:
R:=-10..4.4:
res:=dsolve({ode3,n(0)=1,D(n)(0)=0},numeric,range=R);
plots:-odeplot(res,[n(t),diff(n(t),t)],R,refine=1);

The syntax for Ut(x,0) = 0 is D[2](u)(x,0)=0:

 PDE := diff(u(x,t),x,x)=diff(u(x,t),t,t);
IBC := {u(0,t)=0,u(10,t)=0,u(x,0)=10-x,D[2](u)(x,0) = 0};
 pds := pdsolve(PDE,IBC,numeric);

#And with the piecewise defined initial condition for u:

IBC := {u(0,t)=0,u(10,t)=0,u(x,0)=piecewise(x<5,x,10-x),D[2](u)(x,0) = 0};
pds := pdsolve(PDE,IBC,numeric,spacestep=.01);
pds:-animate(t=5);

You must have an old version of Maple. Many versions ago you had to do as follows (and that obviously will still work):

shift:=(f::procedure)->subs(ff=f,(x->ff(x+1))):

With some Maple version between Maple 4 and 11  lexically scoped variables (I think they are called) were introduced.

You have an obvious division by zero problem with your differential equations for X[b,1,N,D] and X[b,1,B,H] (and more).
I removed that problem in a rather crude way: Replacing the zero intial values by epsilon = 0.001.
Also I saw no point in having S[b,1,O] in the equations at all so I took care of that.
Let sys be your system as given by you above.
sys2:=subs((S[b,1,O](t)=2) = NULL,sys): #Removing S[b,1,O](t)=2
sys3:=(eval(sys2,S[b, 1, O](t) = 2)): #Using that S[b,1,O](t)=2
#res:=dsolve(sys3,numeric); #No success
#res(.1);
odes:=select(has,sys3,diff): #The odes
nops(%);
ics:=select(hastype,sys3,function(0)); #The initial conditions
nops(ics);
for k in odes do k end do; #Taking a closer look
eps:=0.001: #Substitute for zero
#ics2:=lhs~(ics)=~eps; #No need to do this for 11 of the functions
ics3:=remove(has,ics,X[b,1,B,H]);#The problematic one is X[b,1,B,H]
res1:=dsolve(odes union ics3 union {X[b,1,B,H](0)=eps},numeric,maxfun=10^5);
#res1:=dsolve(odes union ics2,numeric,maxfun=10^5);
res1(1); #Testing
##Now you can plot either one of the graphs using odeplot or do them all at once.
##They have widely different values eventually so plotting all in one plot may not be a good idea.
##However, here it is:
depvar:=convert(indets(rhs~(odes),function(identical(t))),list);
plots:-odeplot(res1,[seq([t,depvar[i]],i=1..12)],0..100,legend=[seq(depvar[i],i=1..12)]);
plots:-odeplot(res1,[seq([t,log10(depvar[i])],i=1..12)],0..100,legend=[seq(depvar[i],i=1..12)]);

##############
Less crudely (perhaps) can the problem be overcome by introducing piecewise functions in the 4 differentiall equations where the division by zero occurs.
odes1,odes2:=selectremove(has,odes,{diff(X[b,1,ND](t),t),diff(X[b,1,B,H](t),t),diff(S[b,1,ND](t),t),diff(S[b,1,S](t),t)});
odes1a:=subs(X[b, 1, S](t)/X[b, 1, B, H](t)=piecewise(t<eps,1,X[b, 1, S](t)/X[b, 1, B, H](t)),odes1);
#Now using the zero initial conditions:
res2:=dsolve(odes1a union odes2 union ics,numeric,maxfun=10^5);
#Also you could use the following instead:
piecewise(X[b,1,B,H](t)<eps,1,X[b, 1, S](t)/X[b, 1, B, H](t))

The problem is that f(k); returns an error, when k is just a name.
So either use unevaluation quotes or use the procedural version:
restart;
f:=k->fsolve(x^2+k*x-1=0,x)[1];
f(1);
plot('f(k)',k=0..1);
plot(f,0..1);

You can go to the help page for type and then to
?type,array
A := array(1..2, 1..2, [[1, 3], [1/2, 5]]);
a:=array(1..4,[4,3,2,1]);
type(a,'array'(1,numeric));
type(A,'array'(2,numeric));
type(a,'array'(1,integer));

The result you get by pdsolve(E); is obtained by separation of variables, so can in principle be used to get a series solution as in familiar textbook examples.
By using
pdsolve(E,build);
you will see the building blocks. They look frightening.
So numerical solution is the way to go, but requires concrete values for the parameters and initial and boundary conditions.
Example.
#All parameters have been taken to be 1:
res:=pdsolve(eval(E,{EI=1,L=1,S0=1,p=1,Fv=1}),{x(0,t)=0,x(1,t)=sin(t),x(y,0)=y*(1-y),D[2](x)(y,0)=0},numeric);
res:-plot3d(t=0..0.1);

Here is a not exactly perfect, but yet workable solution:

restart;
FS:=Sum(A[n]*sin(n*Pi*x/L),n=1..infinity); #Capital Pi
u:=combine(sin(m*Pi*x/L)*FS);
applyop(int,1,u,x=-L..L) assuming m::posint ;
simplify(%); #Obviously ignoring the special case n=m
Now consider the special case:
op(1,u);
int(eval(op(1,u),n=m),x=-L..L) assuming m::posint;
#Now taking your function f(x)=x and doing the same operation of multiplying followed by integration:
int(x*sin(m*Pi*x/L),x=-L..L) assuming m::posint;
cf:=eval(solve(%%=%,{A[m]}),m=n);
fs:=subs(cf,FS); #The Fourier series
plot(eval(fs,{infinity=20,L=1}),x=-1..1); #Seems that we got it right!
## You may also try:
value(fs);
plot(eval(%,L=1),x=-1..1);

T:=fsolve((XX(t)-XX(0))^2+(YY(t)-YY(0))^2,t=1..26);
XX(T)-XX(0);
YY(T)-YY(0);
plot([XX,YY,0..T]);
plots:-odeplot(ds,[X(t),Y(t)],0..T,numpoints=500);
plots:-odeplot(ds,[X(t),Y(t)],0..T,frames=20,numpoints=500);

##########
Alternatively you could use events:
##Version 1:
ds := dsolve(odse, numeric, maxfun = 0, abserr = .1^10, relerr = .1^10, minstep = .1^10,events=[[[(X(t)-xf)^2+(Y(t)-yf)^2-1e-6,t>1],halt]]);
ds(30);
##Version 2:
ds := dsolve({odse[],d(0)=0}, numeric, maxfun = 0, abserr = .1^10, relerr = .1^10, minstep = .1^10,discrete_variables=[d(t)::float],events=[[[(X(t)-xf)^2+(Y(t)-yf)^2-1e-6,t>1],d(t)=t]]);
ds(30);

You could try this:

Q:=piecewise(x(t)<=z(t),8,10);
eq1 := diff(x(t), t) = 3*x(t)-1;
eq2 := diff(y(t), t) = y(t)+Q; #Q not Q(t)
eq3 := diff(z(t), t) = z(t);
eqs := {eq1, eq2, eq3};

dsolve(eqs); #Notice the form of the output: A list of 3 sets!
#Compare with
dsolve({eq1,eq3}); #One set!
#Continue with (as an example)
xz_sol:=dsolve({eq1,eq3} union {x(0)=1,z(0)=2});
eval(sol[3],xz_sol);
dsolve(% union {y(0)=0});
y_sol:=combine(expand(value(%)));
plot(eval(Q,xz_sol),t=0..1);

###  Numerical solution is easier:

res:=dsolve(eqs union {x(0)=1,z(0)=2,y(0)=0},numeric);
plots:-odeplot(res,[[t,x(t)],[t,y(t)],[t,z(t)]],0..1);
plots:-odeplot(res,[t,Q],0..1);

To get started go to the help page in Maple 12:
?dsolve,numeric
Take a  look at the first example given at the bottom. It is linear, but that is irrelevant. You can change it as you like.