NDSolveValue::ntdvdae and NDSolveValue::icfail while solving system of PDEs

Question

I am trying to solve a system of PDEs, but I get the error NDSolveValue::ntdvdae:

Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations

After that Mathematica runs for a while and then spits out another error NDSolveValue::icfail:

Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.

Here is the code:

l0 = -1;
r0 = 1;
G = 6.67 10^-11;
c = 1;
c0 = 299792458;
P[t_, r_] := \[Epsilon][t, r]/3;
ic = {R[0, r] == r, Derivative[1, 0][R][0, r] == c, 
   Derivative[1, 0][\[Lambda]][0, r] == 
    1, \[Lambda][0, r] == (-(1/9) r + 10/9) l0, \[Nu][0, 
     r] == (-(1/9) r + 10/9) ((4 Pi)/3 10)^(-1/4), \[Epsilon][0, 
     r] == (4 Pi)/3 10};
bc = {R[t, r0] == r0 + c t, 
   R[t, 10] == 10 + c t, \[Lambda][t, r0] == l0, \[Lambda][t, 10] == 
    0, \[Nu][t, r0] == ((4 Pi)/3 10)^(-1/4), \[Nu][t, 10] == 
    0, \[Epsilon][t, r0] == (4 Pi)/3 10};
eq1 = D[\[Epsilon][t, r], t] R[t, 
      r] + (P[t, r] + \[Epsilon][t, r]) (2 D[R[t, r], t] + 
       R[t, r] D[\[Lambda][t, r], t]) == 0;
eq2 = (P[t, r] + \[Epsilon][t, r]) D[\[Nu][t, r], r] + 
    D[P[t, r], r] == 0;
eq3 = D[\[Nu][t, r], r] D[R[t, r], t] + 
    D[R[t, r], r] D[\[Lambda][t, r], t] - D[D[R[t, r], t], r] == 0;
eq5 = R[t, 
     r] (Exp[2 \[Nu][t, 
          r]] (-((D[\[Lambda][t, r] - \[Nu][t, r], 
              r]) (D[R[t, r], r] + R[t, r] D[\[Nu][t, r], r])) + 
         D[R[t, r], {r, 2}] + R[t, r] D[\[Nu][t, r], {r, 2}]) - 
      Exp[2 \[Lambda][t, 
          r]] (((D[\[Lambda][t, r] - \[Nu][t, r], t]) (D[R[t, r], t] +
              R[t, r] D[\[Lambda][t, r], t])) + D[R[t, r], {t, 2}] + 
         R[t, r] D[\[Lambda][t, r], {t, 2}])) == (8 \[Pi] G)/
    c0^4 Exp[2 \[Lambda][t, r] + 2 \[Nu][t, r]] P[t, r];
{en, Rad, lam, nu} = 
 NDSolveValue[{eq1, eq2, eq3,(*eq4,*)eq5, ic, bc}, {\[Epsilon], 
   R, \[Lambda], \[Nu]}, {t, 0, 10}, {r, r0, 10}]

Answer 1

eq2 can be solved symbolically.

seq2 = DSolve[{eq2, bc[[7]]}, ϵ[t, r], {t, r}][[1, 1]]
(* ϵ[t, r] -> 40/3 E^(4 γ[t, 1] - 4 γ[t, r]) Pi *)

Now, apply the ϵ[0, r] initial condition.

seq2 /. t -> 0 /. (ic[[6]] /. Equal -> Rule);
Solve[% /. Rule -> Equal, γ[0, r]][[1, 1]] /. C[1] -> 0
(* γ[0, r] -> γ[0, 1] *)

which is inconsistent with ic[[5]].

γ[0, r] == ((3/(5 Pi))^(1/4) (10/9 - r/9))/2^(3/4)

It is unlikely that the error messages reported in the question are caused by this inconsistency. Nonetheless, it is not possible to proceed until the inconsistency is resolved.