Here is a (now revised) solution which is very fast also for the rosenbrock method. It takes about the same time as  using lsode[adamsfull] when we use the same solution procedure (solproc) below. On my computer about 0.4s.
##This is a slightly revised version. The first one had (erroneously) X[1](t) instead of just X[1] in the solution procedure. That prevented the use of hardware floats, apparently.
It seems essential for the performance of both methods when used with solproc that solve is not used to isolate the derivatives. The isolation is pretty trivial for this system, but solve does other things in the process as is illustrated below.
When not using the procedural version `DEtools/convertsys` is called. That procedure appears to use solve to isolate the derivatives. The resulting rewriting of the system may be one reason for the slow performance of rosenbrock. It is a curious fact that the nonprocedural version runs very fast when using lsode[adamsfull] though.

restart;
with(plots):
myPathName:="C:/Users/Bruger/Downloads":
read cat( myPathName, "/MC.m");
MC2:=subs(V[1]=X[2],MC): #Using X[1] and X[2] as the elements of vector X
sys,ics:=selectremove(has,MC2,diff): #ics is the set of initial conditions
lhs~(sys);
has(rhs~(sys),diff);
######## So the derivatives are essentially isolated on the left.
############# We shall avoid using the following way of isolating the derivatives.
############# It makes adamsfull go real slow.
sys2:=solve(sys,{diff(X[1](t),t),diff(X[2](t),t)}):
nops~(rhs~(sys2));
nops~(rhs~(sys));
whattype~(rhs~(sys));
#############
odeX,odeVtemp:=op(op~([selectremove(has,sys,diff(X[1](t),t))])):
odeX;
odeV:=diff(X[2](t),t)=rhs(odeVtemp)/.2345666667: ## Instead of using solve
############# The solution procedure:
R1,R2:=op(subs({X[1](t)=X[1],X[2](t)=X[2]},rhs~([odeX,odeV]))): #### Changed from first version
solproc:=subs(RHS1=R1,RHS2=R2,(proc(n,t,X,XP) XP[1]:=RHS1; XP[2]:=RHS2 end proc)): #### Changed from first version
ini:=subs(ics,array([X[1](0),X[2](0)]));
#infolevel[`dsolve/numeric`]:=5:
#sol:=dsolve(number=2,procedure=solproc,start=0,initial=ini,numeric,procvars=[X[1](t),X[2](t)],method=lsode[adamsfull]);
sol:=dsolve(number=2,procedure=solproc,start=0,initial=ini,numeric,procvars=[X[1](t),X[2](t)],method=rosenbrock,maxfun=0);
t0:=time():
odeplot(sol, [[t, X[1](t)],[t, X[2](t)]], 0..9, size=[1800,default]);
time()-t0; #lsode[adamsfull] and rosenbrock about 0.4s
data:=plottools:-getdata(%%);
MX:=data[1][-1];
MV:=data[2][-1];
T:=MX[..,1]: #Times
MXV:=<MX|MV[..,2]>; # t,X,V matrix                                                                       
plot([MX,MV],size=[1800,default]);
N:=numelems(T);
GoodRows:=zip((u,v)-> if verify(u, 7.8..8.2, 'interval') then v end if, convert(T,list), [seq(1..N)]);
printf("\n What odeplot gives\n");
printf("------------------------------------------------\n");
printf("        t               X(t)           V(t)\n");
printf("------------------------------------------------\n");

for k in GoodRows do
   printf("%-15.12f  %-15.12f  %-15.12f\n", seq(MXV[k,n], n=1..3))
end do;
print():
####RESULT (rosenbrock):
 What odeplot gives
------------------------------------------------
        t               X(t)           V(t)
------------------------------------------------
7.830000000000   0.006594922004   0.002815875171
7.875000000000   0.006724509403   0.002960841451
7.920000000000   0.006862302569   0.003162937657
7.965000000000   0.007009462363   0.003381227259
8.010000000000   0.007165994274   0.003556203210
8.055000000000   0.007331727393   0.003835083752
8.100000000000   0.007511890805   0.004281471856
8.145000000000   0.007703617825   0.003952030218
8.190000000000   0.007881895449   0.003971385741

#########
#Now evaluate "pointwise" sol(t) for the times retained

print();
printf("\n What sol(t) gives for the same values of t\n");
printf("------------------------------------------------\n");
printf("        t               X(t)           V(t)\n");
printf("------------------------------------------------\n");

for MyTime in T[GoodRows] do
   MySol := map(u -> rhs(u), sol(MyTime)):
   printf("%-15.12f  %-15.12f  %-15.12f\n", seq(MySol[n], n=1..3))
end do:
print():
#### RESULT (rosenbrock)
 What sol(t) gives for the same values of t
------------------------------------------------
        t               X(t)           V(t)
------------------------------------------------
7.830000000000   0.006594922004   0.002815875171
7.875000000000   0.006724509403   0.002960841451
7.920000000000   0.006862302569   0.003162937657
7.965000000000   0.007009462363   0.003381227259
8.010000000000   0.007165994274   0.003556203210
8.055000000000   0.007331727393   0.003835083752
8.100000000000   0.007511890805   0.004281471856
8.145000000000   0.007703617825   0.003952030218
8.190000000000   0.007881895449   0.003971385741

MaplePrimes16-04-03odesys_procII.mw
## I posted a question about the (undesirable in this case) expansion used by solve in solving a system of equations:
http://mapleprimes.com/questions/211105-Undesired-Expansion-In-Solving?sq=211105

odeplot plots the result, if you make the correction mentioned in my comment to your question

restart;
sys := {((D@@2)(phi))(t) = -2*(D(phi))(t)*(D(theta))(t)*cos(theta(t))/sin(theta(t)),
((D@@2)(theta))(t) = (D(phi))(t)^2*cos(theta(t))*sin(theta(t))-9.8*sin(theta(t))};
ics:={theta(0) = (1/2)*Pi, (D(theta))(0) = 0, phi(0) = (1/2)*Pi, (D(phi))(0) = 1};
res:=dsolve(sys union ics,numeric);
X:= sin(theta(t))*cos(phi(t)):
Y:= sin(theta(t))*sin(phi(t)):
Z:= cos(theta(t)):
plots:-odeplot(res,[X,Y,Z],0..10,labels=[x,y,z]);
## An animation:
plots:-odeplot(res,[X,Y,Z],0..10,frames=100,labels=[x,y,z]);

restart;
integral := Int(x^n*ln(x)^n, x);
IntegrationTools:-Change(integral,x=exp(t));
simplify(%) assuming n>=0,t::real;
value(%);
res:=subs(t=ln(x),%);
diff(res,x) - IntegrationTools:-GetIntegrand(integral);
simplify(%);

Here is an example:

restart;
K:=[seq(LinearAlgebra:-RandomVector(3),i=1..5)];
sort(K,(u,v)->evalb(u[3]<v[3]));

You have a boundary value problem for an ode, not a pde.
Notice that Re means the real part in Maple. I changed it to RE.
You have two conflicting boundary conditions: D(f)(-1)=0, D(f)(-1)=1. I changed the latter to f(-1)=1.

restart;
eq1:=diff(f(eta),eta,eta,eta,eta)+2*(epsilon/(1+epsilon*theta(eta)))^2*diff(f(eta),eta,eta)*((diff(theta(eta),eta))^2)-(epsilon/(1+epsilon*theta(eta)))*(diff(theta(eta),eta,eta))*(diff(f(eta),eta,eta))=0;
eq2:=diff(theta(eta),eta,eta)+Pr*RE*f(eta)*diff(theta(eta),eta)=0;
RE:=1:Pr:=1:epsilon:=0.25:
###bc:=f(1)=0,D(f)(-1)=0,D(f)(1)=1,D(f)(-1)=1,theta(-1)=0,theta(1)=0;#
bc:=f(1)=0,D(f)(-1)=0,D(f)(1)=1,f(-1)=1,theta(-1)=0,theta(1)=0;
##
res:=dsolve({eq1,eq2,bc},numeric);
plots:-odeplot(res,[eta,f(eta)],-1..1);

Besides the problem Carl points out, be aware that implicitplot3d plots the points for which (in this case) ex1(x,t,z) = 0. But ex1(x,t,z) is never zero. As is well known the exp function is never zero even for complex input.
A syntactical aside:
The syntax would be
implicitplot3d( ex1(x,t,z), x = -10..0, t = 0..10, z = 0..10, etc,etc);

For the second question you can use foldl:

foldl(f,a,b,c);
A:=LinearAlgebra:-RandomMatrix(2,2);
B:=LinearAlgebra:-RandomMatrix(3,2);
C:=LinearAlgebra:-RandomMatrix(1,2);
foldl(LinearAlgebra:-KroneckerProduct,A,B,C);
#This gives the same as
LinearAlgebra:-KroneckerProduct(A,B);
LinearAlgebra:-KroneckerProduct(%,C);
##########################

For the first question take a look at
?LinearAlgebra:-IntersectionBasis

Maybe something like this:
u:=sin(x(t))*cos(y(t))*diff(x(t),t)*diff(y(t),t);
indets(u,name);
w:=evalindets(u,function(identical~({x(t),y(t)})),s->op(0,s)(freeze(op(s))));
var:=indets(w,name) minus {t};
mtaylor(w,var,4);
thaw(%);

You are giving 2 initial conditions, but your ode is only of first order. Moreover these two conditions don't agree with your ode1:
eval(ode1,{diff(z(t),t)=0,z(t)=0}); # result 814.8208660 = 0.
Anyway, you can only use z(0)=0.

First solve for diff(z(t),t): There are two solutions. Pick the one you want:

ode2:=solve(ode1,{diff(z(t),t)});
#On my Maple the rising balloon corresponds to the first:
res:=dsolve(ode2[1] union {z(0) = 0},numeric);
plots:-odeplot(res,[t,z(t)],0..5000);

Please notice that gamma in Maple is Euler's constant. In recent Maple versions you can start with local gamma to prevent this fact from bothering you. Otherwise change the name.
I start with a list version of your system.
restart;
local gamma;
sys:=[2*R^2*(diff(w(t), t))*Pi*Omega*h*rho+R^2*(diff(u(t), t, t))*Pi*h*rho-R^2*u(t)*Pi*Omega^2*h*rho = 0, R^2*(diff(v(t), t, t))*Pi*h*rho = 0, -2*R^2*(diff(u(t), t))*Pi*Omega*h*rho+R^2*(diff(w(t), t, t))*Pi*h*rho-R^2*w(t)*Pi*Omega^2*h*rho = 0, (1/4)*R^4*Pi*(diff(alpha(t), t, t))*h*rho+(1/12)*R^2*Pi*(diff(alpha(t), t, t))*h^3*rho+(1/6)*R^2*Pi*(diff(gamma(t), t))*Omega*h^3*rho-(1/12)*R^2*Pi*alpha(t)*Omega^2*h^3*rho = 0, (1/2)*R^4*Pi*(diff(beta(t), t, t))*h*rho-(1/2)*R^4*Pi*beta(t)*Omega^2*h*rho = 0, (1/4)*R^4*Pi*(diff(gamma(t), t, t))*h*rho+(1/12)*R^2*Pi*(diff(gamma(t), t, t))*h^3*rho-(1/6)*R^2*Pi*(diff(alpha(t), t))*Omega*h^3*rho-(1/12)*R^2*Pi*gamma(t)*Omega^2*h^3*rho = 0];
sysD:=convert(sys,D);
idts:=indets(sys,function(identical(t)));
X:=convert(idts,list);
Xn:=map2(op,0,X);
X1:=D~(Xn)(t);
X2:=(D@@2)~(Xn)(t);
M,r:=LinearAlgebra:-GenerateMatrix(sysD,X2);
C,r1:=LinearAlgebra:-GenerateMatrix(convert(-r,list),X1);
K,F:=LinearAlgebra:-GenerateMatrix(convert(-r1,list),X);
M.<X2>+C.<X1>+K.<X> =~F;
### Coming to think about it, you only need one call to GenerateMatrix:
Y:=op~([X2,X1,X]);
A,rA:=LinearAlgebra:-GenerateMatrix(sysD,Y);
A[..,1..6]; #This is M
M; #Check with first method
A[..,7..12]; # This is C
C;
A[..,13..18]; # This is K
K;
rA; #This is F

In Maple input (=1D) you cannot have a<b<c, you must split: a<b, b<c in assuming:
solve(eq,gamma0(t)) assuming -Pi/2 < gamma[1](t) - theta[1](t) - psi[1](t),gamma[1](t) - theta[1](t) - psi[1](t)<Pi/2;

This seems to work:
p:=proc(u) local T;
  if type(u,`+`) then T:=[op(u)] else T:=[u] end if;
  evalindets[flat](T,`*`,op)
end proc;
u:=2*t*exp(t) + 2*t;
p(u); #Works OK
##But the following doesn't because of the product inside sin
p(sin(2*t)); #Error
## A revised version freezes functions:
q:=proc(u) local T,F;
  if type(u,`+`) then T:=[op(u)] else T:=[u] end if;
  F:=evalindets(T,function,freeze);
  thaw(evalindets[flat](F,`*`,op))
end proc;
q(u); #OK
q(sin(2*t)); OK

Assuming that your 2-D input can be interpreted to mean the same as this corrected 1-D version

ode := mu__1*E__1(x)-mu__2*(u(x)-d__1*E__1(x))/d__2+diff(gamma1*E__1(x)-gamma2*(u(x)-d__1*E__1(x))/d__2, x) = 0;

isc := E__1(0) = 0;

and assuming that u(x) is to be considered a known function, then the syntax is:

res:=dsolve({isc, ode},E__1(x));

You will notice that the result res contains an inert integral (Int, grey in print). Using value on res won't help since Maple cannot do the integral unless u is known.
Suppose u=sin, i.e. u(x) = sin(x) for all x, then
eval(res,u=sin);
value(%);
# will give you a finished result. But you can get it on a nicer form:
combine(expand(%));
resf:=collect(%,exp,factor);

tau*f(t) is not an operand in expr[2]. I think it is as simple as that.
My mistake: Apparently not quite. See bottom.
Try the similar example with I replaced by 3:
conv1 := (x) -> eval(x, tau*f(t) = f(-t)):     # Using eval
conv2 := (x) -> algsubs(tau*f(t) = f(-t),x):   # Using algsubs
expr := Vector([
     tau*f(t),
   3*tau*f(t)
]);
conv1(expr),
conv2(expr);
##
Instead you can do
conv1 := (x) -> eval(x, f(t) = f(-t)/tau):
#############
The help for eval:
Usually, the left-hand sides of the equations will be variable names being evaluated.  However, if the left-hand side x of an equation is not a simple name, a syntactic substitution may be attempted which has the same limitations as with subs.
Compare:
eval(a*b*c,a*b=4);
eval(a*b*3,a*b=4);
eval(a*b*I,a*b=4);
Similarly:
subs(a*b=4,a*b*c);
subs(a*b=4,a*b*3);
subs(a*b=4,a*b*I);

J := Int((x^2+2*x+1+(3*x+1)*sqrt(x+ln(x)))/(x*sqrt(x+ln(x))*(x+sqrt(x+ln(x)))), x);
solve(sqrt(x+ln(x))=u,x);
IntegrationTools:-Change(J,x=%);
simplify(%) assuming u>0;
value(%);
res:=subs(u=sqrt(x+ln(x)),%); #Result 2*ln(LambertW(exp(x+ln(x)))+sqrt(x+ln(x)))+2*sqrt(x+ln(x))
##Note added: #After I saw Axel's result I tried simplify and got his simpler version:
res:=simplify(subs(u=sqrt(x+ln(x)),%)) assuming x>0;
## Checking the result:
diff(res,x);
simplify(%) assuming x>0;
IntegrationTools:-GetIntegrand(J);
simplify(%%-%); #0

Note: I made the assumption u > 0. That means that x > LambertW(1) = .5671432904 to 10 digits.
solve(x+ln(x)=0,x); #LambertW(1)