# Step size is effectively zero / Stiff system by NDsolve when the parameter in system of ODEs is varying

As a very new user of Mathematica, i keep encountering errors as shown below when i try to solve my systems of ode with phi1 being value around or greater than 4

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

    (*Version of not transforming a,b,phi*)
Clear["Global*"]

(*constants*)
phi1 = 5; (*error occurs when phi1 is large*)
KcKbRatio = 2;
DbDaRatio = 2;
DcDaRatio = 1;
Cref = 1;
lambda = 2;
CrefLambdaRatio = 1/2;
aEnd = 0.95;
bEnd = 0.8;
cEnd = 0.6;
alpha = 1 - cEnd*CrefLambdaRatio -
KcKbRatio*CrefLambdaRatio/(KcKbRatio + 1)/DcDaRatio*aEnd;
del = $MachineEpsilon; (*sidenotes:when lambda>Cref,graphs tend to look normal*) (*the system of ode and bcs*) ode1 = a''[x] + 2/x*a'[x] - phi1^2*(1 + KcKbRatio)*Phi[x]*a[x] == 0; ode2 = b''[x] + 2/x*b'[x] + phi1^2/DbDaRatio*Phi[x]*a[x] == 0; ode3 = Phi''[x] + 2/x*Phi'[x] - KcKbRatio*phi1^2*CrefLambdaRatio/DcDaRatio*Phi[x]*a[x] == 0; ode = {ode1, ode2, ode3}; bcs = {a'[del] == 0, a[1] == aEnd, b'[del] == 0, b[1] == bEnd, Phi'[del] == 0, Phi[1] == 1 - cEnd*CrefLambdaRatio}; (*ndsolve*) {afunc, bfunc, Phifunc} = NDSolveValue[{ode, bcs}, {a, b, Phi}, {x, del, 1}]; (*Plot*) Plot[{afunc[x], bfunc[x], Phifunc[x], KcKbRatio*Cref/(KcKbRatio + 1)/DcDaRatio/lambda*afunc[x] + alpha(*, Evaluate[{a[x],b[x]}/.exactsol]*)}, {x, del, 1}, PlotLegends -> "Expressions"]  Previously, I asked a similar question with same ode system in this post, but to no avail as well. I had also tried transforming the solution term a, b, Phi in my code by using transformation such as a_transformed = Sqrt[a] and a_transformed = Log[a], but it failed to work as well. Below is my code with transformation. (*Version of trying to transform a,b,phi*) Clear["Global*"] (*constants*) phi1 = 12; KcKbRatio = 2; DbDaRatio = 2; DcDaRatio = 1; Cref = 1; lambda = 2; CrefLambdaRatio = 1/2; aEnd = 0.95; bEnd = 0.8; cEnd = 0.6; alpha = 1 - cEnd*CrefLambdaRatio - KcKbRatio*CrefLambdaRatio/(KcKbRatio + 1)/DcDaRatio*aEnd; del =$MachineEpsilon;
(*sidenotes:when lambda>Cref,graphs tend to look normal*)
(*the system of ode and bcs*)
ode1 = a''[x] + 2/x*a'[x] - phi1^2*(1 + KcKbRatio)*Phi[x]*a[x] == 0;
ode2 = b''[x] + 2/x*b'[x] + phi1^2/DbDaRatio*Phi[x]*a[x] == 0;
ode3 = Phi''[x] + 2/x*Phi'[x] -
KcKbRatio*phi1^2*CrefLambdaRatio/DcDaRatio*Phi[x]*a[x] == 0;

ode = {ode1 /. {a -> (Exp[A[#]] &), Phi -> (Exp[PHI[#]] &)},
ode2 /. {a -> (Exp[A[#]] &), Phi -> (Exp[PHI[#]] &),
b -> (Exp[B[#]] &)},
ode3 /. {a -> (Exp[A[#]] &), Phi -> (Exp[PHI[#]] &)}};
bcs = {A'[del] == 0, A[1] == Log[aEnd], B'[del] == 0,
B[1] == Log[bEnd], PHI'[del] == 0,
PHI[1] == Log[1 - cEnd*CrefLambdaRatio]};
(*ndsolve*)
{Afunc, Bfunc, PHIfunc} =
NDSolveValue[{ode, bcs}, {A, B, PHI}, {x, del, 1}];

(*Plot*)
Plot[{Exp[Afunc[x]], Exp[Bfunc[x]], Exp[(PHIfunc[x])],
KcKbRatio*Cref/(KcKbRatio + 1)/DcDaRatio/lambda*Exp[Afunc[x]] +
alpha(*,Evaluate[{a[x],b[x]}/.exactsol]*)}, {x, del, 1},
PlotLegends -> "Expressions"]

• I run your code (the first one you show above) on V 13.01 and it works with no errors and no warnings? !Mathematica graphics Which code are you asking about? The first one or the second one? Commented Mar 20, 2022 at 16:43
• I get NDSolveValue::ndsz only with the second code and NDSolveValue[..] returns the input: When those two things happen, it's because the shooting method failed. If there is a solution to the BVP, you have to give the method good "StartingInitialConditions" (see this tutorial). If you have some insight from the problem in which this BVP arises, you might have an idea what a reasonable guess for the "StartingInitialConditions" might be. Commented Mar 21, 2022 at 3:20
• @Nasser The first code will have error if phi1 is large, eg: phi1> 4. Second code is just the transformation of ODEs of first code. Commented Mar 21, 2022 at 5:31
• @MichaelE2 As i want to determine the effect of phi1 on this system of ODEs, the range of suitable StartingInitialConditions would vary for different values of phi1. May I ask is it possible to automate the process of solving without encountering this error? Commented Mar 21, 2022 at 5:38
• Try this: Get rid of the derivative BCs and use Method -> "FiniteElement" in the NDSolve` call. -- Does the result look correct? Commented Mar 21, 2022 at 5:50