Plotting the InterpolatingFunction from NDSolve [closed]

I want to plot the solution of a differential equation which i solve it numerically with NDSolve. Here’s the code:

A := 2.23818*^8;
Zh := 1*^11;
a := 1.496*^11;
M := 1.33*^20;
NDSolve[{(1 + (4 A t^2)/(r[t]^2 (r[t]^2 – Zh^2)) +
Zh^2/(r[t]^2 – Zh^2) + (4 A Zh^4 t^2)/(
r[t]^2 (r[t]^2 – Zh^2)^3) + (8 A Zh^2 t^2)/(
r[t]^2 (r[t]^2 – Zh^2)^2))*
r'[t]^2 + ((-4 A t)/(r[t] (r[t]^2 – Zh^2)) – (4 A Zh^2 t)/(
r[t] (r[t]^2 – Zh^2)^2)) r'[t] + A/(r[t]^2 – Zh^2) – (2 M)/
r[t] + M/a == 0, r[0] == 1.496*^11}, r, {t, 0, 3.15*^5}]

And mathematica answers:

{{r->InterpolatingFunction[{{0.,`, 315000.`}},<>]},{r->InterpolatingFunction[{{0.,`, 315000.`}},<>]}

I want to plot the answer as a function of “t”. what is the code?

=================

=================

2 Answers
2

=================

Here is one way to do it.

With[{A = 2.23818*^8, Zh = 1*^11, a = 1.496*^11, M = 1.33*^20},
{rLo, rHi} =
NDSolve[
{(1 +
(4 A t^2)/(r[t]^2 (r[t]^2 – Zh^2)) + Zh^2/(r[t]^2 – Zh^2) +
(4 A Zh^4 t^2)/(r[t]^2 (r[t]^2 – Zh^2)^3) +
(8 A Zh^2 t^2)/(r[t]^2 (r[t]^2 – Zh^2)^2)) r'[t]^2 +
((-4 A t)/(r[t] (r[t]^2 – Zh^2)) – (4 A Zh^2 t) /
(r[t] (r[t]^2 – Zh^2)^2)) r'[t] +
A/(r[t]^2 – Zh^2) – (2 M)/r[t] + M/a == 0,
r[0] == 1.496*^11},
r, {t, 0, 3.15*^5}][[All, 1, 2]]];

Plot[{rHi[t], rLo[t]}, {t, 0, 3.15*^5}]

Update: If you use r[t] instead of r as the second argument of NDSolve you get the desired functions directly:

sol2 = NDSolve[{(1 + (4 A t^2)/(r[t]^2 (r[t]^2 – Zh^2)) + Zh^2/(r[t]^2 – Zh^2) +
(4 A Zh^4 t^2)/(r[t]^2 (r[t]^2 – Zh^2)^3) +
(8 A Zh^2 t^2)/(r[t]^2 (r[t]^2 – Zh^2)^2))* r'[t]^2 +
((-4 A t)/(r[t] (r[t]^2 – Zh^2)) –
(4 A Zh^2 t)/(r[t] (r[t]^2 – Zh^2)^2)) r'[t] +
A/(r[t]^2 – Zh^2) – (2 M)/r[t] + M/a == 0, r[0] == 1.496*^11},
r[t], {t, 0, 3.15*^5}];

funcs2 = r[t] /. sol2

{InterpolatingFunction[{{0., 315000.}}, <>][t], InterpolatingFunction[{{0., 315000.}}, <>][t]}

So you can use funcs2 directly when plotting:

Plot[funcs2, {t, 0, 3.15*^5}]

Original post:

ClearAll[r, t, A, Zh, a, M]
A = 2.23818*^8;
Zh = 1*^11;
a = 1.496*^11;
M = 1.33*^20;
sol = NDSolve[{(1 + (4 A t^2)/(r[t]^2 (r[t]^2 – Zh^2)) + Zh^2/(r[t]^2 – Zh^2) +
(4 A Zh^4 t^2)/(r[t]^2 (r[t]^2 – Zh^2)^3) +
(8 A Zh^2 t^2)/(r[t]^2 (r[t]^2 – Zh^2)^2))* r'[t]^2 +
((-4 A t)/(r[t] (r[t]^2 – Zh^2)) –
(4 A Zh^2 t)/(r[t] (r[t]^2 – Zh^2)^2)) r'[t] +
A/(r[t]^2 – Zh^2) – (2 M)/r[t] + M/a == 0, r[0] == 1.496*^11},
r, {t, 0, 3.15*^5}];

funcs = r /. sol;

Plot[Evaluate[Through@funcs@t], {t, 0, 3.15*^5}]

How it works:

funcs is a list of two pure functions:

funcs

{InterpolatingFunction[{{0., 315000.}}, <>], InterpolatingFunction[{{0., 315000.}},<>]}

funcs[t]

{InterpolatingFunction[{{0., 315000.}}, <>], InterpolatingFunction[{{0., 315000.}},<>]}[t]

We need to push the argument t inside the {..} to to get the functions we need to plot. This can be done in a number of ways:

r[t] /. sol

or

{funcs[[1]][t], funcs[[2]][t]}

or

#[t] & /@ funcs

or

Through[funcs[t]]

all give

{InterpolatingFunction[{{0., 315000.}}, <>][t], InterpolatingFunction[{{0., 315000.}}, <>][t]}