# Extracting data from NDSolve result

Normally I wouldn't ask such question before even read the material from the Wolfram Documentation about what I'm trying to do but this time I have a pretty tigh dead end. I apologize for that. So, what I need to do is to extract a data from a NDSolve integration to use it in different ways. I will show my whole code and make a precise description for what I need. The code is the following:

AU = 1.49597871*10^13;
n = 0.5*10^15;
cg = 6.67384/10^8;
kb = 1.3806488/10^16;
mh = 1.0078250321*(1.660538921/10^24);
mf = 1.98892*10^33;
Te[ϕ_] := 2200 + 440*Cos[ϕ];
Ra[ϕ_] := 1.9 + 0.1*Cos[ϕ + Pi];
t[r_, ϕ_] := Te[ϕ]*(1 - Sqrt[1 - (Ra[ϕ]/r)^2])^(1/4);

k1[r_, ϕ_] := (6.99/10^14)*(t[r, ϕ]/300)^2.8*Exp[-1950/t[r, ϕ]];
k2[r_, ϕ_] := (1.59/10^11)*(t[r, ϕ]/300)^1.2*Exp[-9610/t[r, ϕ]];
k17[r_, ϕ_] := (3.14/10^13)*(t[r, ϕ]/300)^2.7*Exp[-3150/t[r, ϕ]];
k18[r_, ϕ_] := (2.05/10^12)*(t[r, ϕ]/300)^1.52*Exp[-1736/t[r, ϕ]];
k62[r_, ϕ_] := (1.77/10^11)*Exp[178/t[r, ϕ]];
k63[r_, ϕ_] := (1.85/10^11)*(t[r, ϕ]/300)^0.95*Exp[-8571/t[r, ϕ]];
k94[r_, ϕ_] := (1.65/10^12)*(t[r, ϕ]/300)^1.14*Exp[-50/t[r, ϕ]];
k138[r_, ϕ_] := (5.94/10^17)*(t[r, ϕ]/300)^0.17*Exp[65.9/t[r, ϕ]];
k141[r_, ϕ_] := (1/10^11)*Exp[-4800/t[r, ϕ]];
k143[r_, ϕ_] := (1/10^11)*Exp[-4800/t[r, ϕ]];
k144[r_, ϕ_] := 1/10^9;

r1 = k1[r, ϕ]*nH[r]*nOH[r];
r2 = k2[r, ϕ]*nH[r]*nH2O[r];
r17 = k17[r, ϕ]*nH2[r]*nO[r];
r18 = k18[r, ϕ]*nH2[r]*nOH[r];
r62 = k62[r, ϕ]*nO[r]*nOH[r];
r63 = k63[r, ϕ]*nO[r]*nH2O[r];
r94 = k94[r, ϕ]*nOH[r]*nOH[r];
r138 = k138[r, ϕ]*nTi[r]*nO[r];
r141 = k141[r, ϕ]*nTi[r]*nOH[r];
r143 = k143[r, ϕ]*nTi[r]*nH2O[r];
r144 = k144[r, ϕ]*nTiO[r]*nH2O[r];

dens[r_, ϕ_] := n*Exp[(-((2*cg*mf*mh)/(kb*Te[ϕ]*Ra[ϕ]*AU)))*
(1 - (Ra[ϕ]/r)^(1/2))];
dldens[r_, ϕ_] := D[dens[r, ϕ], r]/dens[r, ϕ];

eqns = {
(dldens[r, ϕ]*nH[r] + Derivative[1][nH][r])/dens[r, ϕ] ==
-r1 + r17 + r18 - r2 + r62,
(dldens[r, ϕ]*nOH[r] + Derivative[1][nOH][r])/dens[r, ϕ] ==
-r1 - r141 + r17 - r18 + r2 - r62 + r63 - r94,
(dldens[r, ϕ]*nO[r] + Derivative[1][nO][r])/dens[r, ϕ] ==
r1 - r138 - r62 - r63 + r94,
(dldens[r, ϕ]*nH2[r] + Derivative[1][nH2][r])/dens[r, ϕ] ==
r1 - r17 - r18 + r2,
(dldens[r, ϕ]*nH2O[r] + Derivative[1][nH2O][r])/dens[r, ϕ] ==
-r143 - r144 + r18 - r2 + r94,
(dldens[r, ϕ]*nTi[r] + Derivative[1][nTi][r])/dens[r, ϕ] ==
-r138 - r141 - r143 - r144,
(dldens[r, ϕ]*nTiO[r] + Derivative[1][nTiO][r])/dens[r, ϕ] ==
r138 + r141 + r143 - r144};

ic = {nH[2] == 2/10^7, nOH[2] == 9.8/10^7, nO[2] == 2/10^5, nH2[2] == 5/10,
nH2O[2] == 8/10^6, nTi[2] == 1.96/10^7, nTiO[2] == 2/10^8};

sol = Table[
NDSolveValue[{eqns, ic} /. ϕ -> ϕi, {nH, nOH, nO, nH2, nH2O, nTi, nTiO},
{r, 2, 6}, MaxStepSize -> 0.001, MaxSteps -> Infinity,
AccuracyGoal -> 0, PrecisionGoal -> 0, Method -> {"BDF"}],
{ϕi, 0, 2*Pi, Pi/5}];


Question: how could I extract the data from the solution to a especific value of r and ϕ?

Let's say, for example, that I want to know nOH[6] for all values of ϕ, so

nOH[6] for ϕ=0
nOH[6] for ϕ=Pi/5
nOH[6] for ϕ=Pi


and so on. Is there some way to do it?

Moreover, is there some way to export the data from the solution where I can obtain an arrangement like the following?

Thank you very much.

• I get several NDSolveValue::ndnum errors when I execute your code. Perhaps you forgot to post something? Dec 21, 2015 at 20:46
• maybe it should be cg*mf*mh ..? Dec 21, 2015 at 21:16
• I think this is what you want: TableForm[ Table[{r, sol[[1, 2]][r]} , {r, {2, 2.1, 2.2, 2.3}}] , TableHeadings -> {None, {r, y}}] Dec 21, 2015 at 21:29
• @garej yes, it look reasonable! And yes, cgmfmh= cg * mf * mh. I'm sorry, didn't see it got the messy when I wrote it here. Dec 21, 2015 at 21:30
• I get sol[[1, 2]][6] -> 258767. Andre, you should edit the question to fix that cgmfmh typo. Dec 21, 2015 at 21:33

In your case you get the matrix of solutions. You can extract them by parts like sol[[2,1]][6] but it is not convenient since you have to think which number stands for what combination of variables: sol[[phi, element]][r].

It is comfortable to set up a Manipulate for the task. There are of course many ways to do that.

Manipulate[
sol;
expr = Thread[List[Range[up - 1, up, 0.1], Table[sol[[phi, element]][r],
{r, Range[up - 1, up, 0.1]}]]];

export = Grid[{
(* for Export first row is better to make as {"\[Phi]=", phi/.phirange} *)
{Item["\[CapitalPhi]=" Dynamic[phi /. phirange],
Alignment -> {Left, Center}, Frame -> {{True, False}}], Nothing},
{"r", element /. list},
##
}, Frame -> All, ItemSize -> {{10, 10}, {3, {1.6}}}] & @@ expr,

{{up, 3, "Range r"}, 3, 6, 0.1, Appearance -> "Labeled"},
Delimiter,
{{element, 2}, list, SetterBar},
Delimiter,
{{phi, 1}, phirange, SetterBar},

Delimiter,
(* change Print[export] for Export["file__.xls",export]*)
Button["Export", Print[export]],

ControlPlacement -> Left,

Initialization :> (
list =
Thread[Range[7] -> {"nH", "nOH", "nO", "nH2", "nH2O", "nTi", "nTiO"}];
phirange =
Thread[Range[Length@Range[0, 2 Pi, Pi/5]] -> Range[0, 2 Pi, Pi/5]];
)
]


Result:

Addenum: With solution as sol some, I would recommend:

Manipulate[
sol;(*Dimensions[sol] gives {11,7}; elements of sol are function heads*)
expr =
TableForm[
Table[sol[[phi, element]][r],
{r, Range[2, up, 0.1]},
{phi, Range[Length@Range[0, p, Pi/5]]}],
TableHeadings -> {Range[2, 6, 0.1], Range[0, 2 Pi, Pi/5]}
];
Framed[expr]
,
Button["Print", Print[expr]],
(* you may change Print[expr] for Export["file.ext", expr]*)

{{element, 1}, {1 -> "nH", 2 -> "nOH", 3 -> "nO", 4 -> "nH2",
5 -> "nH2O", 6 -> "nTi", 7 -> "nTiO"}},
{{up, 3, "Range r"}, 3, 6, 0.1},
{{p, Pi/5, "Range phi"}, Pi/5, 2 Pi, Pi/5}]