How to increase efficiency of Mathematica code?

Question

I have a code that runs as it is supposed to, but it takes around 8 hours to reach the solution with the step size necessary to find agreement between the analytical and computational solutions. As I add complexity to the code, this time is expected to increase.

I'm looking for suggestions on how to increase the efficiency of the code in order to decrease the computational time for the nH2comp interpolating function in the code below, which represents the computational solution in the r and z directions.

So far I have considered:

• Parallelizing: the dependent nature of the table I am generating precludes this solution.
• Larger Step Size: the solution doesn't converge everywhere. I am currently running through a couple of tests to see if a happy medium exists; the largest step size I've found producing an acceptable result is rstep=0.1 and zstep=0.01, which takes the aforementioned 8 hours.

The code below is an example code set to rstep=zstep=1, this takes around 30 seconds for the computational portion (ie.nH2comp). There is a graph included where convergence/ non-convergence can be observed.

Thank you for the help!

(*Defining functions and constants*)
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));
Pref = 1;
a = 4;
zH2 = 10^5;

(*Define r and z step sizes*)
rstep = 1.;
zstep = 1.;

(*Computational Solution*)
Clear[nH2]
Timing[nH2comp = Interpolation[Flatten[Table[{{r, zb},
      nH2[zb] /. 
       NDSolve[{nH2'[z] == -zH2*nH2[z]*(P[r, z]/Pref )^a, 
          nH2[-100] == 1}, {nH2}, {z, -100, zb}][[1]]},
     {r, 0.001, 100, rstep}, {zb, -100, 50, zstep}], 1]]]

(*Analytical Solution*)
rstable = Table[rj, {rj, 0, 100, rstep}];
zstable = Table[zj, {zj, -100, 50, zstep}];
IsTable = 
  Table[0, {j, 1, Dimensions[rstable][[1]]}, {k, 1, 
    Dimensions[zstable][[1]]}];
Timing[Do[
  Do[IsTable[[j, k]] = 
    IsTable[[j, k - 1]] + 
     P[rstable[[j]] + 0.0001, zstable[[k]]]^
      a (zstable[[k]] - zstable[[k - 1]])/Pref^a,
   {k, 2, Dimensions[zstable][[1]]}], {j, 1, 
   Dimensions[rstable][[1]]}];
 IsFunTable = 
  Table[{rstable[[j]], zstable[[k]], IsTable[[j, k]]}, {j, 1, 
    Dimensions[rstable][[1]]}, {k, 1, Dimensions[zstable][[1]]}];
 IsFun = Interpolation[Flatten[IsFunTable, 1]];]

nH2ana [r_, z_] := Exp[-zH2 IsFun[r, z]]

(*Plotting *)
r = {5, 2, 1.0, 0.1};
rplot = Labeled[
  Plot[{nH2comp[r[[1]], z], nH2comp[r[[2]], z], nH2comp[r[[3]], z], 
    nH2comp[r[[4]], z], nH2ana[r[[1]], z], nH2ana[r[[2]], z], 
    nH2ana[r[[3]], z], nH2ana[r[[4]], z]}, {z, -8, 8}, 
   PlotLabel -> {{"r positions:", r}, {"Step size:", rstep "=rstep", 
      zstep "=zstep"}}, PlotRange -> {0, 2}, 
   PlotLegends -> {"Numeric" r[[1]], "Numeric" r[[2]], 
     "Numeric" r[[3]], "Numeric" r[[4]], "Analytic" r[[1]], 
     "Analytic" r[[2]], "Analytic" r[[3]], "Analytic" r[[4]]}, 
   Axes -> True, AxesLabel -> {"z", "nH2"}, 
   PlotStyle -> {{Red, Thick}, {Orange, Thick}, {Brown, 
      Thick}, {Purple, Thick}, {Red, Thick, Dashed}, {Orange, Thick, 
      Dashed}, {Brown, Thick, Dashed}, {Purple, Thick, Dashed}}, 
   PlotRange -> All], {{"nH2 for" , zH2 "=zH2"}}, {{Top, Left}}]

Answer 1

It's hard to understand what you're after, but if it is to solve your differential equation over z for different r, it should be much simpler. Your loop ends up running NDSolve over and over again, there's no need to Interpolate the results yourself, and let NDSolve handle the step sizes automatically.

r = 5;
sol = NDSolve[{nH2'[z] == -zH2*nH2[z]*(P[r, z]/Pref)^a, nH2[-100] == 1},
  nH2, {z, -100, 50}][[1]];
Plot[nH2[z] /. sol, {z, -100, 50}]

EDIT: Based on OP's comment and following @SjoerdSmit's advice to use ParametericNDSolve, how's this look?

Clear[r]
sol = ParametricNDSolve[{nH2'[z] == -zH2*nH2[z]*(P[r, z]/Pref)^a, nH2[-100] == 1},
  nH2, {z, -100, 50}, {r}];

denplot = 
 DensityPlot[nH2[r][z] /. sol, {z, -8, 8}, {r, 0, 15}, 
  FrameLabel -> {"z", "r"}, ImageSize -> Medium, 
  ColorFunction -> "SunsetColors", 
  PlotLegends -> BarLegend[{"SunsetColors", {0, 1}}]]