It's because `NDSolve` is using different method for these 2 problems. For the first problem, `Shooting` method, which is for nonlinear boundary value problem of ordinaray differential equation, is used; while for the second problem, `FiniteElement` method, which supports nonlinear problem as of _v12_, is used. Different solvers of `NDSolve` are not equally powerful. For the second problem, we can obtain a result by adjusting the option a little. (The result is not perfect though): sol = First@ NDSolve[{eqns, ic}, {r[R, t], TRR[R, t], γRR}, {R, 0.1, 1}, {t, 0, 1}, InitialSeeding -> {r[R,t] == 3, TRR[R,t] == 3, γRR[R,t] == 2}, Method -> {FiniteElement, MeshOptions -> MaxCellMeasure -> 0.05}] PlotSol@sol [![enter image description here][1]][1] "Then, is it possible to avoid using pure `FiniteElement`?" Sadly the answer seems to be no, at least now, because adding `Method -> {MethodOfLines, TemporalVariable -> t}` causes `tvic` warning, which seems to indicate `NDSolve` is having difficulty to notice the system can actually be discretized to a DAE system of `t`. Still, if we discretize the system outside of `NDSolve` (I'll use [`pdetoode`](http://mathematica.stackexchange.com/a/127997/1871) for the task): γθθ[R_, t_] := 4 bleft = 1/10; bright = 1; With[{r = r[R, t], TRR = TRR[R, t], γRR = γRR[R, t], γθθ = γθθ[R, t]}, eqns = {r D[r, R] == γRR γθθ R, r^4 R γθθ D[TRR, R] == 2 γRR (r^4 - R^4 γθθ^4), D[γRR, t] == 0}; bc = {r == 0.1 /. R -> bleft, TRR == 0 /. R -> bright}; ic = γRR == 3 /. t -> 0;]; Block[{γRR}, γRR[R_, t_] := 3; solinitlst = NDSolveValue[{eqns // Most, bc}, {r[R, t], TRR[R, t]}, {R, 0.1, 1}]]; icadd = {r[R, 0], TRR[R, 0]} == solinitlst // Thread; tend = 1; points = 100; grid = Array[# &, points, {bleft, bright}]; difforder = 2; (* Definition of pdetoode isn't included in this post, please find it in the link above. *) ptoofunc = pdetoode[{r, TRR, γRR}[R, t], t, grid, difforder]; ode = {D[{Rest@ptoofunc@eqns[[1]], Most@ptoofunc@eqns[[2]]}, t], ptoofunc@eqns[[3]]}; odebc = bc // ptoofunc; odeic = {ic // ptoofunc, icadd[[1]] // ptoofunc // Rest, icadd[[2]] // ptoofunc // Most}; var = Outer[#@#2 &, {r, TRR, γRR}, grid]; sollst = NDSolveValue[{ode, odebc, odeic}, var, {t, 0, tend}]; {solr, solTRR, solγRR} = rebuild[#, grid, 2] & /@ sollst; Plot[{solr, solTRR, solγRR}[R, tend] // Through // Evaluate, {R, bleft, bright}] `NDSolve` will manage to find an accurate enough result: [![enter image description here][2]][2] The approach is excessive for this toy example, but can be useful for more complicated problem. #Remark 1. Without the `D[…, t]` in definition of `ode`, the program will still find a solution, but it's not that accurate, maybe because `NDSolve` can handle differential equation involving derivative of `t` better. 2. `solinitlst` is necessary, probably because the DAE solver can't find accurate enough initial conditions automatically in this case. [1]: https://i.sstatic.net/6felw.png [2]: https://i.sstatic.net/PVJo5.png