Following this question, I’m trying to use NDSolve to solve a system of algebraic and ordinary differential equations:

$ \frac{\partial P}{\partial x}=f\frac{\rho \nu^2}{2D_t} \tag{1}$

$ \frac{\nu^2}{2}+k_3P^{\left(1-\frac{1}{\gamma}\right)}=k_1\tag{2}$

for detailed background you may see this and this.

mu = 1.568*10^-5 rho = 1.225 r = 0.01 P1 = 6*10^5 Lt = 0.01 R0 = 286.9 gamma = 1.4 cP = 1004 T0 = 300 rho1 = P1/(T0*R0) k2 = P1/rho1^gamma k3 = cP*k2^(1/gamma)/R0 k1 = k3*P1^(1 - 1/gamma) DT[c_] := 2*Sqrt[r^2 - (r - c)^2] RE[v_, c_] := rho*v*DT[c]/mu FRHP[v_, c_] := 64/RE[v, c] FRBP[v_, c_] := 0.3164/RE[v, c]^0.25 FRKN[v_, c_] :=  NSolve[(1/Sqrt[f] == (2*Log[RE[v, c]*Sqrt[f]] - 0.8)), f][[1, 1,  2]] // Quiet FR[v_, c_] := If[RE[v, c] < 2320, FRHP[v, c], If[RE[v, c] < 10^5, FRBP[v, c], FRKN]] sol[c_] := NDSolve[{D[P[x], x] == FR[V[x], c]*rho*V[x]^2/(2*DT[c]),       V[x]^2/2 + k3*P[x]^(1 - 1/gamma) == k1, P[0] == P1}, {V, P}, {x,       0, Lt}][[1, 1, 2]] // Quiet 

Mathematica finishes the computations without any errors or warnings but when I plot the result it is empty:

Plot3D[Evaluate[{V[x], P[x]} /. First[sol[c]]], {x, 0, Lt}, {c, 0, r}] 

I expect to see two plots P and V versus x and c. I would appreciate if you could help me know where I’m making a mistake and how to solve it.