I have four coupled nonlinear partial differential equations with one spatial coordinate “mu” and one temporal “t”. These equations are linked to other algebraic equations and together are:

```
G[mu_, t_] = 4 \[Pi]*rho[mu, t]*r[mu, t]^2*D[r[mu, t], mu];(*MW39*)
w[mu_, t_] = 1 + ep[mu, t]/c^2 + p[mu, t]/(rho[mu, t]*c^2);(*MW41*)
a[mu_, t_] = 1/w[mu, t];
ep[mu_, t_] = k*rho[mu, t]^(\[Gamma] - 1)/(\[Gamma] - 1);
p[mu_, t_] = (\[Gamma] - 1) ep[mu, t]*rho[mu, t];(*MW40*)
equt[mu_,
t_] = -a[mu,
t] (4 \[Pi]*r[mu, t]^2*G[mu, t]/w[mu, t]*D[p[mu, t], mu] + (
m[mu, t]*gr)/
r[mu, t]^2 + (4 \[Pi]*gr)/c^2*p[mu, t]*r[mu, t]);(*MW33*)
eqrt[mu_, t_] = a[mu, t]*u[mu, t];(*MW34*)
eqmm[mu_, t_] =
4 \[Pi]*rho[mu, t]*(1 + ep[mu, t]/c^2)*
r[mu, t]^2 D[r[mu, t], mu];(*MW38*)
eqrhort[mu_, t_] = -a[mu, t]*rho[mu, t]*r[mu, t]^2 D[u[mu, t], mu]/
D[r[mu, t],
mu];(*MW35-new variable rho*r^2->Subscript[(rho*r^2), t]=eqrhot*)
\
(*equations and conditions*)
eqs = {D[u[mu, t], t] == equt[mu, t], D[r[mu, t], t] == eqrt[mu, t],
D[m[mu, t], mu] == eqmm[mu, t],
D[rho[mu, t]*r[mu, t]^2, t] == eqrhort[mu, t]};
bcon = {DirichletCondition[u[mu, t] == 0., mu == dmu],
DirichletCondition[r[mu, t] == r0[dmu], mu == dmu],
DirichletCondition[m[mu, t] == fm0[dmu], mu == dmu],
DirichletCondition[rho[mu, t] == frho0[dmu], mu == dmu]};
incon = {u[mu, 0] == 0., r[mu, 0] == r0[mu], m[mu, 0] == fm0[mu],
rho[mu, 0] == frho0[mu]};
(*solution*)
NDSolveValue[{eqs, incon, bcon}, {u, r, m, rho}, {mu, dmu, mumax}, {t,
0, 0.1}]
```

Initial functions r0[mu], fm0[mu] and frho0[mu] are interpolated functions coming from numerical solution of stationary problem. Result of this solution are error messages:

```
NDSolveValue::femcnsd: The PDE coefficient -((6.674*10^-8 m[mu])/r[mu]^2)-1.15712*10^-17 r[mu] rho[mu]^(5/3)-3.26355*10^23 r[mu]^4 rho[mu]^(2/3) (r^\[Prime])[mu] (rho^\[Prime])[mu] does not evaluate to a numeric scalar at the coordinate {2.08798*10^34}; it evaluated to Indeterminate instead.
NDSolveValue::femcnsd: The PDE coefficient -((6.674*10^-8 m[mu])/r[mu]^2)-1.15712*10^-17 r[mu] rho[mu]^(5/3)-3.26355*10^23 r[mu]^4 rho[mu]^(2/3) (r^\[Prime])[mu] (rho^\[Prime])[mu] does not evaluate to a numeric scalar at the coordinate {2.08798*10^34}; it evaluated to Indeterminate instead.
```

Unfortunately, I don’t know where is the problem (whole concept, method or…). The problem appears ever in half value of endpoint of integration (mumax/2), doesn’t matter what “mumax” is. I’m able to draw (and evaluate in all point of the range) all defined functions in the initial time without problems.

Thank you for help or suggestions.

PS: I’m new here if something is misspelt, marked or unlisted. Please notify me.

Thank you.