Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Here is an example of the system that I am running,

eqns = {y'[t] == y[t] (x[t] - 1), x'[t] == x[t] (a - y[t]), x[0] == 1,y[0] == 2.7};
tmax = 100;
pfun = ParametricNDSolveValue[eqns, Integrate[y[t], {t, 0, tmax}] , {t, 0, tmax},{a}];
Plot[pfun[a], {a, 0, 1}, AxesLabel -> {"a", "y"}]

and for a particular range of the parameter, "a", I get numerical giberish as shown here:

enter image description here

To check what is going on in this region I looked at the timeseries of x and y for different values of "a". If you plot example timeseries for "a" values higher than the numerical noise you get nice oscillating curves enter image description here

This was found using this code:

tmax = 100;
amin = 1;
amax = 3;
aint = 1;
eqns1[a_] := {y'[t]==y[t] (x[t] - 1), x'[t]==x[t] (a - y[t]), x[0]==1, y[0]==2.7};
pn1 = ParametricNDSolveValue[eqns1[a], x, {t, 0, tmax}, {a}]
pn2 = ParametricNDSolveValue[eqns1[a], y, {t, 0, tmax}, {a}]

plt1 =Plot[Evaluate[Table[pn1[a][t],{a,amin,amax,aint}]], {t,0,tmax},PlotRange-> All];
plt2 =Plot[Evaluate[Table[pn2[a][t],{a,amin,amax,aint}]], {t,0,tmax},PlotLegends->  
Range[amin, amax, aint], PlotRange -> All];
GraphicsRow[{plt1, plt2}]

However, to get the timeseries to plot for the a values around the noise, i.e. between 0 and 0.1, I had to add Workingprecision-> 20 and Method-> "LSODA" into pn1 and pn2 (into the ParametricNDSolveValue options). Which nicely renders the curves enter image description here

So the problem is that these options are not necessary for all the "a" values when I plot pfun above. In fact, if I include these options into my real system (which is larger) the whole plot takes exponentially longer to finish. So I would like to turn on these options for only some "a" values and not others, when I make a plot of pfun. Or use another trick you might know.

I would greatly appreciate any suggestions on how to do this, or other comments. Thank you.

share|improve this question
up vote 1 down vote accepted

You can switch between the options by defining a function that calls one form of ParametricNDSolve if a is in a certain range, and otherwise calls the other form.

eqns = {y'[t] == y[t] (x[t] - 1), x'[t] == x[t] (a - y[t]), x[0] == 1,
    y[0] == 2.7`20};
tmax = 100;

pfun = ParametricNDSolveValue[eqns, Integrate[y[t], {t, 0, tmax}], {t, 0, tmax}, {a}];
pfun2 = ParametricNDSolveValue[eqns, Integrate[y[t], {t, 0, tmax}], {t, 0, tmax}, {a}, 
   WorkingPrecision -> 20, Method -> "LSODA"];

myp[a_] /; 0.04 < a < 0.1 := pfun2[a];
myp[a_] := pfun[a];

You still may have some work to do on your DE system, however:

Plot[myp[a], {a, 0, 1}, AxesLabel -> {"a", "y"}]

Mathematica graphics

For example:


ParametricNDSolveValue::ndsz: At t$159261 == 83.68333550162991842139594420675849673339`20., step size is effectively zero; singularity or stiff system suspected. >>

InterpolatingFunction::dmvali: The integration endpoint 100 in dimension 1 lies outside the range of data in the interpolating function. Extrapolation will be used. >>

(* Out= 0.*10^63 *)
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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