# Plotting NDSolve results when you change one of the parameters

I have the following ODE that I have solved with NDSolve:

y'[x] == -l*a*(y[x]^2 - c*y[x] - p*x^3*Exp[-2*x])/x^2


where $l = 1.508*10^{13}$, $a = 5*10^{-9}$, $c = 10^{-11}$, and $p = 1.038*10^{-5}$ are constants using the code

sol1 = NDSolve[{y'[x] == -l*a*(y[x]^2 - c*y[x] - p*x^3*Exp[-2*x])/x^2,
y == 10^(-5)}, y, {x, 10, 100}, MaxStepSize -> 0.001].


But now that I've done this, I want to vary $a$. Specifically, I want to graph $y$ as a function of $a$, where $a$ varies from $10^{-7}$ to $10^{-11}$.

I've tried looking at past threads but I haven't seen anyone do what I want to do.

I've tried making sol1 a function of $a$ but that doesn't work because then I've fed NDSolve a non-numerical value.

I've tried solving the ODE as a PDE (not an ideal solution...) by taking $a$ as a variable with the following code

sol2 =
NDSolve[{D[y[x, a], x] == -l*a*(y[x, a]^2 - c*y[x, a] - p*x^3*Exp[-2*x])/x^2,
y[10, a] == 10^(-5)}, y, {x, 10, 100}, {a, 10^-11, 10^-7}, MaxStepSize -> 0.001].


However, I run into an error anyway: "At x == 10., step size is effectively zero; singularity or stiff system suspected."

Does anyone know how to graph the end result of sol1, $y$, as a function of $a$?

• If NDSolve runs fast in your case, you can solve the equation multiple times for differene $a$, retrieve the value of $y(100)$ for each and plot a ListPlot. – corey979 Sep 7 '16 at 18:15
• I could, but it'd really be ideal to get a continuous function. I'm trying to replicate the results of a paper, and their $f$ v $a$ plot is continuous. – Ecclesiastic Sep 7 '16 at 18:21
• But you know that in a computer nothing is contiuous? Even when do Plot[x^2,{x,0,1}] you in fact obtain a set of connected points, not a mathematically continuous graph. You can plot yours for e.g. 100 or 1000 values of $a$. – corey979 Sep 7 '16 at 18:25
• I suppose. But I think the difference would be that if I want to obtain $y$ for, say, $a = 1.5*10^{-11}$ and I only have calculated $y$ for the points $a = 10^{-11}, 2*10^{-11}, ..., 10^{-7}$ beforehand, then I will have to go calculate $y$ again, whereas if I want the value of $x^2$ for any $x$, I can obtain it much more easily. So that's the sort of sense in which I'd like the function to be continuous. It seems like that should be possible. – Ecclesiastic Sep 7 '16 at 18:33
• ParametricNDSolve often is used for this purpose. – bbgodfrey Sep 7 '16 at 19:07

As mentioned by bbgodfrey in his comment, you can use ParametricNDSolve, or otherwise ParametricNDSolveValue:

l = 1.508 10^13;
c = 10^(-11);
p = 1.038 10^(-5);

eqs = {y'[x] == -l*a*(y[x]^2 - c*y[x] - p*x^3*Exp[-2*x])/x^2, y == 10^(-5)};

sol = ParametricNDSolveValue[eqs, y, {x, 10, 100}, {a}] Plot[sol[a], {a, 10^-11, 10^-7}]
` • Thanks so much! This is what I needed. As a quick note, I actually copied the l value wrong--it's $3.01593*10^{21}$. This results in bad behavior for the above graph but this can be fixed by adding AccuracyGoal -> 20 in the ParametricNDSolveValue. – Ecclesiastic Sep 8 '16 at 5:19