# Pulling Results from NDSolve within a Table

I am still learning the intricacies of Mathematica. I have an NDSolve within which I am trying to vary a specific parameter (I understand that ParameticNDSolve is likely a more efficient tool--the results it is producing are physically incorrect, I am hoping NDSolve will allow for the correct solution). To do this, I have a Table on the outside of the NDSolve. The result is a table with interpolating functions. I need help figuring out how to extract these results from the table so that I can graph them. Thank you for the help!

Note: Once I can extract values from this table I plan to build two more tables so as to vary additional variables (d and zH below).


P[r_, z_] := (2. r ((2 - (4 r)/((1 + r)^2 +
z^2)) EllipticK[(4 r)/((1 + r)^2 + z^2)] -
2 EllipticE[(4 r)/((1 + r)^2 + z^2)]))/(\[Pi] ((4 r)/((1 +
r)^2 + z^2)) \[Sqrt]((1 + r)^2 + z^2));
Zmaxreg = 15;
Zminreg = -Zmaxreg;
Rmaxreg = 54;
Rminreg = 0;

d = 0;
zH = 0;
pref = 1;
a = 4;

Clear[H2, H]
Quiet[
Table[sol = NDSolve[
{(*EQUATIONS *)
D[H2[r, z], z] == -(zH2 + d)*H2[r, z]*(P[r, z]/pref )^a,
D[H[r, z], z] ==
2*d*H2[r, z]*(P[r, z]/pref )^a - zH*H[r, z]*(P[r, z]/pref )^a,

(*INITIAL CONDITIONS*)
H2[r, Zminreg] == 1,
H[r, Zminreg] == 0
},
(*VARIABLES TO SOLVE FOR*)
{H2, H},
(*RANGE TO SOLVE OVER*)
{r, Rminreg, Rmaxreg}, {z, Zminreg, Zmaxreg}
],
{zH2, 1, 10^3, 10}
]
]

H2[r_, z_] = H2[r, z] /. sol[[1]];
H[r_, z_] = H[r, z] /. sol[[1]];

~~~~~

• Please post a minimal working example. Jun 30, 2022 at 0:38
• Well, you already know there's a Part ([[]]), don't you? That's exactly what you need. Please check the document of it more carefully. Jun 30, 2022 at 2:50
• @xzczd I read through the documentation prior to posting and wasn't able to find the answer I needed. Are you able to elaborate on how exactly I would be able to pull out a single interpolating function from the table?
– Bill
Jul 7, 2022 at 22:53
• It's confusing to me you're having difficulty with [[]]. (Especially after seeing you claimed you know how to do that in MATLAB in comment below. It's similar in Mathematica. ) How do you understand the sol[[1]] in your code? You simply copy it from somewhere without understanding it? Anyway, observe the output of the following code: solutionlist = Table[NDSolve[{x'[t]==-y[t]-x[t]^2,y'[t]==2 x[t]-y[t]^3,x[0]==k, y[0]==1}, {x, y}, {t, 0, 1}], {k, 1, 7}]; solutionlist[[3]], then observe solutionlist[[3, 1]], then solutionlist[[3, 1, 1]], then solutionlist[[3, 1, 1, -1]]. Jul 8, 2022 at 0:13

\$Version

(* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *)

Clear["Global*"]

P[r_, z_] := (2 r ((2 - (4 r)/((1 + r)^2 + z^2)) *
EllipticK[(4 r)/((1 + r)^2 + z^2)] -
2 EllipticE[(4 r)/((1 + r)^2 + z^2)]))/
(π ((4 r)/((1 + r)^2 + z^2)) √((1 + r)^2 + z^2));
Zmaxreg = 15;
Zminreg = -Zmaxreg;
Rmaxreg = 54;
Rminreg = 0;

d = 0;
zH = 0;
pref = 1;
a = 4;

sys = {(*EQUATIONS*)
D[H2[r, z], z] == -(zH2 + d)*H2[r, z]*(P[r, z]/pref)^a,
D[H[r, z], z] ==
2*d*H2[r, z]*(P[r, z]/pref)^a - zH*H[r, z]*(P[r, z]/pref)^a,
(*INITIAL CONDITIONS*)
H2[r, Zminreg] == 1, H[r, Zminreg] == 0};


Using NDSolve

Clear[H2, H]
Quiet[
tab = Table[{zH2,
NDSolve[sys,(*VARIABLES TO SOLVE FOR*)
{H2, H},
(*RANGE TO SOLVE OVER*)
{r, Rminreg, Rmaxreg}, {z, Zminreg, Zmaxreg}]},
{zH2, 1, 10^3, 10}]];

Manipulate[
Plot3D[Evaluate[{H2[r, z], H[r, z]} /. tab[[n, 2]]],
{r, Rminreg, Rmaxreg}, {z, Zminreg, Zmaxreg},
AxesLabel -> (Style[#, 14] & /@ {r, z}),
PlotLabel -> Style[StringForm["zH2 = ", tab[[n, 1]]], 14],
PlotLegends -> (Style[#, 14] & /@ {H2, H})],
{{n, 61}, 1, Length[tab], 1, Appearance -> "Labeled"}]


sol = ParametricNDSolve[sys,
(*VARIABLES TO SOLVE FOR*)
{H2, H},
(*RANGE TO SOLVE OVER*)
{r, Rminreg, Rmaxreg}, {z, Zminreg, Zmaxreg},
(*Parameter*)
{zH2}];

Manipulate[
Plot3D[Evaluate[{H2[zH2][r, z], H[zH2][r, z]} /. sol],
{r, Rminreg, Rmaxreg}, {z, Zminreg, Zmaxreg},
AxesLabel -> (Style[#, 14] & /@ {r, z}),
PlotLegends -> (Style[#, 14] & /@ {H2, H})],
{{zH2, 601}, 1, 1000, 10, Appearance -> "Labeled"}]
`

• Thank you for the response! That is helpful. What I'm still not sure about is how to access a specific cell within the table instead of all of them at once. For example, say that I would like to pull only H2 for zH2 = 200. In MatLAB, that would be a simple matter of understanding the shape of the matrix that the data is stored in and calling it. I'm struggling to find the equivalent operation in Mathematica.
– Bill
Jul 7, 2022 at 20:25