# Allow Mathematica to solve the differential equations with larger than 16 digits parameters

Consider the system of differential equations

test = ParametricNDSolve[{x'[t] ==
1.1*(1 - x[t])*a*y[t] - x[t]*a*z[t],
y'[t] == -1.1 (1 - x[t])*b*y[t] - c*y[t] + d*Exp[-t/0.02],
z'[t] == -x[t]*b*z[t] - c*z[t] + d*Exp[-t/0.02], y[0.02] == 0,
z[0.02] == 0, x[0.02] == 0.1}, {x, y, z}, {t, 0.02, 0.5}, {a, b, c,
d}]


If I set the parameters a,b,c,d to

av = 1;
bv = 2.75*10^6;
cv = 4*10^7;
dv = 10^16;


then the solutions are obtained. However, setting, for example, d = 10^37, and a = 10^-23, when calculating the solution

solutiontest1 = x[av, bv, cv, dv] /. test


Mathematica signals about an error:

ParametricNDSolve::ndsz: At t == 0.02, step size is effectively zero; singularity or stiff system suspected.

How to resolve this problem?

• @Rebel, usually, we do not accept answers in this site so quickly; this is so that other users can give other (and potentially more refined) answers. A 24-hour waiting period is common. Commented Nov 14, 2019 at 14:45

As the error says, the system is stiff, so you need to use a Method that is appropriate for stiff systems, either "BDF" or "StiffnessSwitching". From experience, I've found that it usually helps to also increase the WorkingPrecision to some high value and use exact numbers, eg 11/10 instead of 1.1 and so on.

For example, the following works:

a = 1;
b = 275/100*10^6;
c = 4*10^7;
d = 3*10^16;
p0 = 2/100;
test = NDSolve[{x'[t] == 11/10*(1 - x[t])*a*y[t] - x[t]*a*z[t],
y'[t] == -(11/10) (1 - x[t])*b*y[t] - c*y[t] + d*Exp[-t/p0],
z'[t] == -x[t]*b*z[t] - c*z[t] + d*Exp[-t/p0], y[p0] == 0,
z[p0] == 0, x[p0] == 1/10}, {x, y, z}, {t, p0, 1/2},
WorkingPrecision -> 25, Method -> "BDF", MaxSteps -> 10^6][[1]];
LogLogPlot[{x[t], y[t], z[t]} /. test // Abs // Evaluate, {t, 0.02,0.5}]


Alternatively, you can also use "StiffnessSwitching", but I found it makes the solution a bit wiggly.

• The desired value of a parameter d is something of an order of 10^37, while of a parameter a is something like 10^-23. With this value, the error persists. Commented Nov 14, 2019 at 15:15
• You did not specify in the question what is the "desired" value of d. Please post all necessary information, otherwise any answer will be incomplete... Commented Nov 14, 2019 at 15:16
• I have already added it. But initially, I mentioned values d > 10^16, which can be any, not necessarily close to 10^16, so the value 3*10^16 was your assumption for which the answer works. Commented Nov 14, 2019 at 15:18