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: <br>
<pre><code>
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);</code>

<pre><code>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;</pre></code>

<pre><code>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];</pre></code>

<pre><code>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, ϕ];</pre></code>


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


<pre><code>ic = {nH[2] == 2*10^-7, nOH[2] == 9.8*10^-7, 
   nO[2] == 2*10^-5, nH2[2] == 5*10^-1, 
   nH2O[2] == 8*10^-6, nTi[2] == 1.96*10^-7, 
   nTiO[2] == 2*10^-8};</pre></code>                      

<pre><code>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}];</pre></code>
 
Now, 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 <br>
nOH[6] for ϕ=0 <br>
nOH[6] for ϕ=Pi/5 <br>
nOH[6] for ϕ=Pi <br>
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 a arrangement like the following?

[![Table][1]][1]


  [1]: https://i.sstatic.net/9B1Z6.png

Thank you very much.