# Have been getting this error "NDSolve::ndsz: At z == 0.8068293549897082, step size is effectively zero; singularity or stiff system suspected."

Tried googling the error "step size is effectively zero; singularity or stiff system suspected". Haven't found any soln. Any help would mean a lot. Posting the full code.

w0 = -0.5; w1 = 0.05; wm = 0.01; \[CapitalOmega]m0 = 0.23; \
\[CapitalOmega]c0 = 0.78; \[Delta] = 0.01; H = 72;

h1[z_] := \[CapitalOmega]m0 ((1 + z)^(3 (1 + wm - \[Delta])) +
3 \[Delta] (1 + z)^(-3 (1 + w0 - w1))*
E^(w1 (1 + z))*(-3 w1)*(1 + z)^(-3 (w0 - w1 - wm + \[Delta]) -
1)*(-(-3 w1)*(1 + z))^(3 (w0 - w1 - wm + \[Delta]) -
1) Gamma[-(3 (w0 - w1 - wm + \[Delta]) - 1) -
1, -(-3 w1) (1 + z)]) + \[CapitalOmega]c0 (1 + z)^(
3 (1 + w0 - w1)) E^(3 w1 (1 + z));
s1 = NDSolve[{E^(K[z]/3) K'[z] == 3/Sqrt[h1[z]], K[.01] == .01},
K[z], {z, 0, 10}];
f1[z_] := Evaluate[K[z] /. s1];

w00 = -1;
h2[z_] := \[CapitalOmega]m0 ((1 + z)^(3 (1 + wm - \[Delta])) +
3 \[Delta] (1 + z)^(-3 (1 + w00 - w1))*
E^(w1 (1 + z))*(-3 w1)*(1 +
z)^(-3 (w00 - w1 - wm + \[Delta]) -
1)*(-(-3 w1)*(1 + z))^(3 (w00 - w1 - wm + \[Delta]) -
1) Gamma[-(3 (w00 - w1 - wm + \[Delta]) - 1) -
1, -(-3 w1) (1 + z)]) + \[CapitalOmega]c0 (1 + z)^(
3 (1 + w00 - w1)) E^(3 w1 (1 + z));
s2 = NDSolve[{E^(K[z]/3) K'[z] == 3/Sqrt[h2[z]], K[.01] == .01},
K[z], {z, 0, 10}];
f2[z_] := Evaluate[K[z] /. s2];

w000 = -1.5;
h3[z_] := \[CapitalOmega]m0 ((1 + z)^(3 (1 + wm - \[Delta])) +
3 \[Delta] (1 + z)^(-3 (1 + w000 - w1))*
E^(w1 (1 + z))*(-3 w1)*(1 +
z)^(-3 (w000 - w1 - wm + \[Delta]) -
1)*(-(-3 w1)*(1 + z))^(3 (w000 - w1 - wm + \[Delta]) -
1) Gamma[-(3 (w000 - w1 - wm + \[Delta]) - 1) -
1, -(-3 w1) (1 + z)]) + \[CapitalOmega]c0 (1 + z)^(
3 (1 + w000 - w1)) E^(3 w1 (1 + z));
s3 = NDSolve[{E^(K[z]/3) K'[z] == 3/Sqrt[h3[z]], K[.01] == .01},
K[z], {z, 0, 10}];
f3[z_] := Evaluate[K[z] /. s3];

Plot[{f1[z], f2[z], f3[z]}, {z, 0, 5}] 

• What do the last two lines have to do with your problem? Why are they included if they don't matter? Nov 7, 2021 at 12:45
• Didn't quite get your question? Whatever lines I have added all are necessary. Please try to run the code in your machine. P.S- Thank you for yesterday's help. It worked at last. Nov 7, 2021 at 12:49
• The last two lines of code are f3[z_] := Evaluate[K[z] /. s3]; and Plot[{f1[z], f2[z], f3[z]}, {z, 0, 5}]. The problem occurs in the 3rd-to-last line with NDSolve. I don't want to Plot a failed NDSolve. It potentially wastes time, distracts from the real problem, causes lots of error messages, all of which are pointless until the NDSolve code is fixed. I sometimes just skip such questions. -- P.S. You're welcome. Nov 7, 2021 at 13:03
• Actually, you should just include the one NDSolve that fails, not all three. Why make it so difficult for other to find your problem? Nov 7, 2021 at 13:05
• Will keep that in mind from next time. In my machine I am getting problem with the 2nd and 3rd NDSolve. First one worked fine and gave the graph. Nov 7, 2021 at 13:09

The integrand in your second integration seems to be stiff what would require extremeliy small steps. Therefore, use the method "StiffnessSwitching" in the second integration like:

s2 = NDSolve[{E^(K[z]/3) K'[z] == 3/Sqrt[h2[z]], K[.01] == .01},
K[z], {z, 0, 10}, Method -> "StiffnessSwitching"]


• Tried adding  Method -> "StiffnessSwitching", the error got way but still I am getting one graph instead of three. Nov 7, 2021 at 12:51
• @DarthVader Don't think that might be because of the same reason as last time? (Complex numbers.) Also you might want to define f1[z_] := Evaluate[K[z] /. First@s1] Nov 7, 2021 at 13:10
• It seems I got the same complex number problem like yesterday. Will try to rectify it my changing the parameters. Nov 7, 2021 at 13:49

First, don't use K as your own symbol. It is used internally by Mathematica.

You can determine whether it's stiffness (often fixable) or a singularity (a feature which does not need fixing, and can't be "fixed"). Plotting the solution or the norm of the solution might suggest a singularity. Trying a non-stiff solver such as an explicit RK method will test for stiffness.

s2 = NDSolve[{E^(k[z]/3) k'[z] == 3/Sqrt[h2[z]], k[.01] == .01},
k, {z, 0, 10}, Method -> "ExplicitRungeKutta"]

(* {{k -> InterpolatingFunction[{{0., 10.}}, "<>"]}} *)


No stiffness error, so it must be a problem in the (Automatic) LSODA method.

Well, duh, there's an obvious singularity when h2[z] == 0, which occurs at the point z == 0.8... in the error. I should have checked this first, but a little slow this morning, I guess.

Plot[h2[z], {z, 0, 1}]


Since it's a pole of order 1/2, it's integrable, but I don't know of a graceful way to get NDSolve to solve it. The method "ExplicitRungeKutta" must accidentally skip over it due to round-off error. In fact, if we set WorkingPrecision -> 16, we get a divide-by-zero error.

• We are allowed to change the parameters w0, w1, wm, $\delta$. Will see if changing them solves the issue. Meanwhile, if you were able to find something do kindly share. Thank you. Nov 7, 2021 at 13:54