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

Question

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}] ```

Answer 1

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

Answer 2

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.