# How to plot the derivative of NDSolve's interpolation function? [duplicate]

i just get the answer (graph) from code below. Right now i want to get the differentiate of d[t]; d''[t]. but i don't know how to get it. This is my code

m = (E^(-(d[t]^2/(a[t]^2 + 1))))/(a[t]^2 + 1)^(5/2)
n = Sqrt[\[Pi]]*a[t]
V0 = -1
W0 = -2
s =
NDSolve[{c'[
t] == (2*
m)*(W0*(a[t]^4*b[t] - 2*a[t]^2*b[t]*d[t]^2 + a[t]^2*b[t] +
a[t]^2*d[t]*c[t] + d[t]*c[t]) - a[t]^2*d[t]*V0),
d'[t] == (a[t]^2*W0*m*(a[t]^2 - 2*d[t]^2 + 1)) + c[t] ,
a'[t] ==
a[t]*(1 + a[t]^2)^-1*d[t]*W0*
m*(-2*a[t]^4 + 2*a[t]^2*d[t]^2 - a[t]^2 + 1) + 2*a[t]*b[t],
b'[t] ==
1/(2 a[t]^4) - n/(2*Sqrt[2 \[Pi]]*a[t]^3) -
V0*m (a[t]^2 + 2*d[t]^2 + 1) - 2*b[t]^2, c[0] == 0.25,
d[0] == -12, a[0] == 1.0, b[0] == 0}, {d, c, a, b}, {t, 0, 120}]

dsol = d[t] /. s
ddsol = D[d[t], t, t] /. s
ParametricPlot[{dsol, ddsol}, {t, 0, 120}, Frame -> True,
PlotRange -> {{-5, 5}, {-2, 2}},
FrameLabel -> {StyleForm["\[Xi] (center-of-mass) ", FontSize -> 14],
StyleForm["F(\[Xi]) (force)", FontSize -> 14]}]


i get the value of d[t] (graph). My question is, how to differentiate that value; d''[t], next?

(i need to plot a graph, x axis is for d, and for y axis is d''. [i need to get graph as picture below)

• You could use ParametricPlot[] if you want to make a phase portrait. Commented Mar 30, 2017 at 23:53
• ParametricPlot [{d[t], d''[t]}, {t, -5, 5}] ? Commented Mar 31, 2017 at 0:16

The structure of the function, you want to plot is $f(t)=(d(t),\,d''(t))$, which is the parametric equation. For this, as suggested by @J.M., you should use ParametricPlot.

xi = Exp[-d[t]^2/(a[t]^2 + 1)]/(a[t]^2 + 1)^(5/2); N1 = Sqrt[Pi]*A^2*a[t];
W0 = 0.2; A = 1; V0 = 5.86;

eq1 = v'[t] == (2*xi)*(W0*(a[t]^4*b[t] - 2*a[t]^2*b[t]*d[t]^2 + a[t]^2*b[t] +
a[t]^2*d[t]*v[t] + d[t]*v[t]) - a[t]^2*d[t]*V0);

eq2 = d'[t] == a[t]^2*W0*xi*(a[t]^2 - 2*d[t]^2 + 1) + v[t];

eq3 = a'[t] == a[t]*(1/(1 + a[t]^2))*d[t]*W0*
xi*(-2*a[t]^4 + 2*a[t]^2*d[t]^2 - a[t]^2 + 1) + 2*a[t]*b[t];

eq4 = b'[t] == (1/(2*a[t]^4)) - N1/(2*Sqrt[2*Pi]*a[t]^3) -
V0*xi (a[t]^2 + 2*d[t]^2 + 1) - 2*b[t]^2;

sol = First@NDSolve[{eq1, eq2, eq3, eq4, v[0] == 0.25, d[0] == -12, a[0] == 1,
b[0] == 0.0}, {v, d, a, b}, {t, -120, 120}];

dsol = d[t] /. sol;

ddsol = D[d[t], t, t] /. sol;

ParametricPlot[{dsol, ddsol}, {t, -120, 120}, PlotRange -> {{-5, 5}, {-2, 2}}]


• i try it, then i try to change the value of V0=-1, W0=-2, and time from 0 to 120. why i don't get a good graph (from the article)? Commented Mar 31, 2017 at 18:47

m = (E^(-(d[t]^2/(a[t]^2 + 1))))/(a[t]^2 + 1)^(5/2)
n = Sqrt[\[Pi]]

{d[t_], c[t_], a[t_], b[t_]} = {d[t], c[t], a[t], b[t]} /. NDSolve[{c'[
t] == (2*
m)*(2*(a[t]^4*b[t] - 2*a[t]^2*b[t]*d[t]^2 + a[t]^2*b[t] +
a[t]^2*d[t]*c[t] + d[t]*c[t]) - a[t]^2*d[t]*-1),
d'[t] == (a[t]^2*2*m*(a[t]^2 - 2*d[t]^2 + 1)) + c[t],
a'[t] ==
a[t]*(1 + a[t]^2)^-1*d[t]*2*
m*(-2*a[t]^4 + 2*a[t]^2*d[t]^2 - a[t]^2 + 1) + 2*a[t]*b[t],
b'[t] ==
1/(2 a[t]^4) - n/(2*Sqrt[2 \[Pi]]*a[t]^3) +
1*m (a[t]^2 + 2*d[t]^2 + 1) - 2*b[t]^2, c[0] == 0.25,
d[0] == -12, a[0] == 1.0, b[0] == 0}, {d, c, a, b}, {t, 0,
120}][[1]]
d''[1]