# Stiff system suspected error under certain parameter values

My code is working fine until I try to increase the value of parameter phi1, error such as

NDSolveValue::ndsz: At x == 0.23372152690117168, step size is effectively zero; singularity or stiff system suspected.

starts to appear. Previously, I fixed this error in my code by employing this trick. However, this does not work anymore when i increase the value of parameter phi1

Clear["Global*"]
(*constants*)
phi1 = 12;
KcKbRatio = 2;
DbDaRatio = 2;
DcDaRatio = 0.5;
Cref = 1;
lambda = 2;
aEnd = 0.95;
bEnd = 0.8;
cEnd = 0.6;
alpha = 1 - cEnd *Cref/lambda -
KcKbRatio*Cref/(KcKbRatio + 1)/DcDaRatio/lambda*aEnd ;
del = $MachineEpsilon; (*sidenotes: when lambda>Cref, graphs tend to look normal *) (*the system of ode and bcs*) ode = {a''[x] + 2/x*a'[x] - phi1^2*(1 + KcKbRatio)*phi[x]*a[x] == 0, b''[x] + 2/x*b'[x] + phi1^2/DbDaRatio*phi[x]*a[x] == 0, phi''[x] + 2/x*phi'[x] - KcKbRatio*phi1^2*Cref/DcDaRatio/lambda*phi[x]*a[x] == 0}; bcs = {a'[del] == 0, a[1] == aEnd , b'[del] == 0, b[1] == bEnd, phi'[del] == 0, phi[1] == (lambda - cEnd *Cref)/lambda}; (*ndsolve*) {asol, bsol, phisol} = NDSolveValue[{ode, bcs}, {a, b, phi}, {x, del, 1}]; (*Dsolve*) exactsol = DSolve[{a''[x] + 2/x*a'[x] - phi1^2*(1 + KcKbRatio)*(1 - cEnd*Cref/lambda)*a[x] == 0, a'[0] == 0, a[1] == aEnd, b''[x] + 2/x*b'[x] + phi1^2/DbDaRatio*(1 - cEnd*Cref/lambda)*a[x] == 0, b'[0] == 0, b[1] == bEnd}, {a[x], b[x]}, x]; (*Plot*) Plot[{asol[x], bsol[x], phisol[x], KcKbRatio*Cref/(KcKbRatio + 1)/DcDaRatio/lambda*asol[x] + alpha, Evaluate[{a[x], b[x]} /. exactsol]}, {x, del, 1}, PlotLegends -> "Expressions"]  • I quickly tried the code you provided on version 12, and I get the stifness error even for phi1 = 0.01. Is it possible that you missed some part of your code? – kcr Mar 11 at 17:47 • I am using version 13.0 – joe Mar 11 at 17:59 • thanks for that! – kcr Mar 11 at 17:59 • @kcr Does the solution code from this work? – joe Mar 11 at 18:03 • yup, the answer from that link works fine on v12 – kcr Mar 11 at 18:04 ## 1 Answer It works when using a piecwise expression for x -> 0. Since a'[0] == 0  the expression 2 a'[x]/x  tends to 1. Therefore substitute for a and equivalent for b and phi (2 Derivative[1][a][x])/x -> Piecewise[{{1, x <$MachineEpsilon}}, (2 Derivative[1][a][x])/x];

Clear["Global*"]
(*constants*)
phi1 = 12;
KcKbRatio = 2;
DbDaRatio = 2;
DcDaRatio = 0.5;
Cref = 1;
lambda = 2;
aEnd = 0.95;
bEnd = 0.8;
cEnd = 0.6;
alpha = 1 - cEnd*Cref/lambda -
KcKbRatio*Cref/(KcKbRatio + 1)/DcDaRatio/lambda*aEnd;
del = \$MachineEpsilon;
(*sidenotes:when lambda>Cref,graphs tend to look normal*)

(*the system of ode and bcs*)
ode = {a''[x] + 2/x*a'[x] - phi1^2*(1 + KcKbRatio)*phi[x]*a[x] == 0,
b''[x] + 2/x*b'[x] + phi1^2/DbDaRatio*phi[x]*a[x] == 0,
phi''[x] + 2/x*phi'[x] -
KcKbRatio*phi1^2*Cref/DcDaRatio/lambda*phi[x]*a[x] == 0};

bcs = {a'[del] == 0, a[1] == aEnd, b'[del] == 0, b[1] == bEnd,
phi'[del] == 0, phi[1] == (lambda - cEnd*Cref)/lambda};

ode2 = {-432*a[x]*phi[x] + Piecewise[
{{1, x < $$MachineEpsilon}}, (2*Derivative[1][a][x])/x] + Derivative[2][a][x] == 0, 72*a[x]*phi[x] + Piecewise[{{1, x < MachineEpsilon}}, (2*Derivative[1][b][x])/x] + Derivative[2][b][x] == 0, -288*a[x]*phi[x] + Piecewise[{{1, x <$$MachineEpsilon}},
(2*Derivative[1][phi][x])/x] + Derivative[2][phi][x] ==
0};


Method "StiffnessSwitching" produced quite large errors for lhs of ode. Method "BDF" performes better, as seen by the second graph.

ndsol = NDSolve[Join[ode2, bcs], {a, b, phi}, {x, 0, 1},
Method -> "BDF"]

Plot[Evaluate[{a[x], b[x], phi[x]} /. ndsol[[1]]], {x, 0, 1},
PlotStyle -> {Red, Green, Blue}, PlotRange -> All,
GridLines -> {Automatic, {19/20, 4/5, 7/10}}]


Plot[Evaluate[ode[[All, 1]] /. ndsol[[1]]], {x, 0, 1},
PlotStyle -> {Red, Green, Blue}, PlotRange -> .02]


May be a more smooth approach to zero than piecwise could improve solutions near zero.

.

• It is not working for my mathematica 13.0.1 for Microsoft Windows (64-bit)` The following error would occur when i copy this into my mathematica:
– joe
Mar 15 at 15:29