# Plotting the derivative of the result of NDsolve [duplicate]

I have the following set of equations:

X[r_] := (-4 a*A[r]/(r^2 ϕ[r] Sqrt[A[r] B[r]])
D[r^2/Sqrt[A[r] B[r]] D[ϕ[r], r], r] + 2 a*m^2)/(6 α^2);

Scalar1[r_] := β*D[r^2*Sqrt[A[r]/B[r]]*D[ψ[r], r], r]- 3*k*α^2*X[r]*(ϕ[r])^2*r^2*Sqrt[B[r]/A[r]];

TensorTT1[r_] := (A[r] (-B[r] + B[r]^2 + r Derivative[B][r]))/(r^2 B[r]^2) + (3 α^2 X[r] - a*m^2) (ϕ[r])^2/2 -2 a*(D[ϕ[r], r])^2/(
2 B[r]) - β*(D[ψ[r], r])^2 A[r]/(2 B[r]);

TensorRR1[r_] := (A[r] - A[r] B[r] + r Derivative[A][r])/(r^2 A[r]) + (3 α^2 X[r] \[Minus] a*m^2) B[r] (ϕ[r])^2/(
2 A[r]) + 2 a*(D[ϕ[r], r])^2/(2 A[r]) - β*(D[ψ[r], r])^2/2;

Tensorθθ1[r_] := (r (-r B[r] Derivative[A][r]^2 - 2 A[r]^2 Derivative[B][r] +
A[r] (-r Derivative[A][r] Derivative[B][r] +
2 B[r] (Derivative[A][r] + r A''[r]))))/(4 A[r]^2 B[
r]^2) + (3 α^2 X[r] - a*m^2) r^2 (ϕ[r])^2/(
2 A[r]) - 2 a*r^2 (D[ϕ[r], r])^2/(2 A[r] B[r]) + β*
r^2*(D[ψ[r], r])^2/(2 B[r]);


And I am solving it using NDSolve:

sol1 = NDSolve[{Scalar1[r] == 0, TensorTT1[r] == 0, TensorRR1[r] == 0,
Tensorθθ1[r] == 0, A == 1, B == 1, ϕ == 0.01, ψ ==
0.01, ϕ' == 0.1, ψ' == 0.001, A' == 0.001}, {A[r], B[r], ϕ[r], ψ[r]}, {r, 0.1, 100}]


I get a nice clean plot when I try:

Plot[Evaluate[A[r] /. sol1], {r, 0.1, 0.2}, PlotRange -> All]


But I am more interested in looking at the plot of the derivative, but doing

Plot[Evaluate[A'[r] /. sol1], {r, 0.1, 50}, PlotRange -> All]


gives me an empty plot.

• To be more specific, your problem will be resolved if you modify the 2nd argument of NDSolve to {A, B, \[Phi], \[Psi]}. – xzczd Oct 19 '17 at 5:56
• Include A'[r] in your list of functions to be solved for, i.e., NDSolve[eqns, {A[r], A'[r], B[r], ϕ[r], ψ[t]}, {r, .1, 100}]. This should provide a better result than numerically differentiating the interpolating function for A[r]. – Carl Woll Oct 19 '17 at 5:58
• @CarlWoll I think it gives the same result as differentiating the interpolating function. Differentiating an interpolating function is done in a smart way, by recomputing the interpolation coefficients, not by (f[x+delta] - f[x]) / delta with a small delta (or similar). I just tried this out with {ff,fd} = NDSolveValue[..., {f,f'}, ...] and ff' == fd gives True`. – Szabolcs Oct 19 '17 at 19:10