I don't believe that there is any way of getting a formula for the general solution, nor is there a way to find a formula for the ivp {eq, y(0) = 283.8}.

But you can easily solve numerically:

res:= dsolve({eq, y(0) = 283.8},numeric);
plots:-odeplot(res,[x,y(x)],-1500..1465);


####
# Note:
dsolve(eq);
returned nothing (NULL) in Maple 17.02 as well as in Maple 2016 on my computer.

Use fnormal.
See the help page ?fnormal.
But here is a simple example (Digits is default 10 on my computer):
restart;
Digits; #10
u:=a*1e-3+b*1e-9+c*1e-15;
fnormal(u);
fnormal(u,11);
fnormal(u,17);

The D you get in the following is the usual D.
restart;
b:=<b1(t),b2(t)>;
f(b); #Not a container (relevant for elementwise ~)
##The only container in the following is b:
dfdb:=Physics[diff]~(f(b),b); #Output a vector with two identical elements! Is that intended?
dfdb[1];
               D(f)(<b1(t), b2(t)>)  ##Output (same as dfdb[2])

If we now define f as
f:=b->b[1]^2+b[2];

then notice that even now the following doesn't evaluate to anything other than the input:
D(f)(b);
               D(f)(<b1(t), b2(t)>)  ##Output
###
If you make f a function of two variables instead of a function of a vector you get more of what you want:
restart;
b:=<b1(t),b2(t)>;
dfdb:=Physics[diff]~(f(b1(t),b2(t)),b);
f:=(b1,b2)->b1^2+b2;
dfdb;

I think the explanation is that the harm has already been done before evalf gets to the number:
restart;
5^1.25; #10 digits with default settings of Digits
'5^1.25'; #Unevaluation quotes don't help: automatic simplification.
5^(5/4); #Still symbolic

So raise Digits:
Digits:=30:
evalf(5^1.25);

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;