Stiff system suspected error under certain parameter values

Question

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

Answer 1

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.

.