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*)
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, 
  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, 
 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}}]
  • $\begingroup$ Your analytical solution nH2ana looks like numerical one since it based on numerical table IsFunTable and interpolating function IsFun. Then it is not clear what do you compare with on the rplot? $\endgroup$ Aug 16, 2022 at 4:54
  • $\begingroup$ Could you edit your question / title to make it more specific? It seems that the real question is "How to plot how the solution of a differential equation depends on a parameter?" That will help others with similar questions find this in the future. $\endgroup$
    – Chris K
    Aug 16, 2022 at 14:50

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

enter image description here

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

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

enter image description here

  • $\begingroup$ Hello--I ultimately need to see the solution over all r. I'm running it for a slice of r as a way to check if the solution is correct. I thought of using the automatic discretization within the NDSolve, but it results in different lengths, which then makes it difficult to store and use later on. The final plot I'm looking for is something like below: 'denplot = DensityPlot[nH2comp[r, z], {z, -8, 8}, {r, 0, 15}, FrameLabel -> {"z", "r"}, ImageSize -> Medium, ColorFunction -> "SunsetColors", PlotLegends -> BarLegend[{"SunsetColors", {0, 1}}]]' $\endgroup$
    – Bill
    Aug 16, 2022 at 3:16
  • 4
    $\begingroup$ @Bill Have you looked at ParametricNDSolveValue? That's specifically for dealing with NDSolve problems with varying model parameters. $\endgroup$ Aug 16, 2022 at 7:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.