How to solve the following GR PDE?

Question

Someone could help me for solving the following GR PDE? (this code do not work yet)

E0 = 1; G = 1; c = 1; 
R = 5; T = 1/3; 

Cte = Sqrt[(E/Pi/R^3)*(E0/(4/T^2 + 1/R^2))]
Phi = Cte*Exp[(-r/2)*R]*Exp[(-t^2/4)*T^2]; 
pressure = (1/2)*(Exp[2*F[r, t]]*D[Phi, t]^2 + D[Phi, r]^2)

eqrr = ((-c^2)*E^(2*F[r, t]) + 
E^(2*H[r, t])*((-E^(2*F[r, t]))*
    r^2*(3*Derivative[0, 1][H][r, t]^2 + 
      2*Derivative[0, 2][H][r, t]) + (c + 
      c*r*Derivative[1, 0][H][r, t])^2))/
 (E^(2*H[r, t])*(c^2*r^2))
eqthetatheta = (1/c^2)*E^(-2*F[r, t] + 2*H[r, t])* r*((-E^(2*F[r, t]))*
 r*(Derivative[0, 1][F][r, t]^2 + 
   Derivative[0, 1][F][r, t]*Derivative[0, 1][H][r, t] + 
   Derivative[0, 1][H][r, t]^2 + 
         Derivative[0, 2][F][r, t] + Derivative[0, 2][H][r, t]) + 
c^2*((-Derivative[1, 0][F][r, t])*(1 + 
      r*Derivative[1, 0][H][r, t]) + 

   Derivative[1, 0][H][r, t]*(2 + r*Derivative[1, 0][H][r, t]) + 
   r*Derivative[2, 0][H][r, t]))

precis = 30; Tsim = 10*T; Rsim = 10*R; 
NDSolve[{eqrr == 8*Pi*G*pressure, eqthetatheta == 0, F[r, -Tsim] == 0, H[r, -Tsim] == 0, 
    Derivative[1, 0][F][0, t] == 0, Derivative[1, 0][H][0, t] == 0, F[Rsim, t] == 0, 
    H[Rsim, t] == 0}, {F, H}, {t, -Tsim, Tsim}, {r, 0, Rsim}, WorkingPrecision -> precis, 
    Method -> {"MethodOfLines", "SpatialDiscretization" -> 
    {"TensorProductGrid", "MaxPoints" -> 151,"MinPoints" -> 151, 
    "DifferenceOrder" -> "Pseudospectral"}, "TemporalVariable" -> t}]

Answer 1

As I noted in comments above, NDSolve has boundary value problems. It needs boundary conditions at t = -Tsim for the first derivatives of F and H with respect to time, and it cannot have spatial boundary conditions involving derivatives of F with respect to r. In the absence of additional information, I modified the arguments of NDSolve as follows.

sol = NDSolve[{eqrr == 8*Pi*G*pressure, eqthetatheta == 0, 
    F[r, -Tsim] == 0, H[r, -Tsim] == 0, Derivative[1, 0][H][r0, t] == 0,
    F[Rsim, t] == 0, H[Rsim, t] == 0, 
    (D[F[r, t], t] /. t -> -Tsim) == 0, (D[H[r, t], t] /. t -> -Tsim) == 0}, 
    {F, H}, {t, -Tsim, Tsim}, {r, r0, Rsim}, Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 151, 
    "MinPoints" -> 151, "TemporalVariable" -> t}]

with r0 = 1/100 to avoid singularities at r = 0, and Tsim = T to reduce run time to a manageable amount. The output is

Plot3D[Evaluate[{F[r, t], H[r, t]} /. sol], {t, -Tsim, -.22}, {r, r0, Rsim}, 
    PlotRange -> All]

Clearly, the computation is failing near r = r0. To make progress, consider carefully the proper boundary conditions, as opposed to those I simply made up.