Extracting data from NDSolve result

Question

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.

Answer 1

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}]

Special thanks to @george2079 comments.