# Error using Nintegrate within NDsolve

I'm trying to numerically solve a system of differential equations which contains functions defined using NIntegrate:

Sa[y_] := (2*π*y)/(1 - Exp[-2 π*y]);
MSa[z_?NumericQ] := 2/√π NIntegrate[Sa[(z/u)^(1/2)]*u^(1/2)*Exp[-u], {u, 0, +∞}];

Sb[y_] := (2^10*π)/3*y^5/(1 + y^2)^2*Exp[-4 y*ArcCot[y]]/(1 - Exp[-2 π*y]);
MSb[z_?NumericQ] :=  2/√π NIntegrate[Sb[(z/u)^(1/2)]*(u^(1/2)*Exp[z])/(Exp[z + u] - 1)*Exp[-u], {u, 0, +∞}];

Yeq[x_?NumericQ] := c3 x^(3/2) Exp[-x];

fion[z_?NumericQ] :=  2^7/3 NIntegrate[η/(1 + η^2)^2*Exp[-4 η*ArcCot[η]]/(1 - Exp[-2 π*η])*1/(Exp[z (1 + 1/η^2)] - 1), {η, 0, +∞}];

odes := {Y'[z] == (-c1*MSa[z])/z^2 (Y[z]^2 - Yeq[4*z/α^2]^2) - (c1*MSb[z])/z^2 Y[z]^2 + c2*z*fion[z] (S[z] + T[z]), S'[z] == (c1*MSb[z])/4 z^2*Y[z]^2 - c2 z (1 + fion[z]) S[z], T'[z] == (3*c1*MSb[z])/(4*z^2)*Y[z]^2 - c2 z (ca + fion[z]) T[z], Y[(5*α^2)/4] == Yeq, S[(5*α^2)/4] == 1/2*Yeq[10 - (5*α^2)/4], T[(5*α^2)/4] == 3/2*Yeq[10 - (5*α^2)/4]}

sol = NDSolve[odes /. {m -> 1000, α -> 137^-1, Subscript[M, pl] -> 2.44*10^18}, {Y, S, T}, {z, (5*137^-2)/4, 10}]


where c1, c2, c3 and ca are fixed constants given by:

c1 := 2.13 Subscript[M, pl] m^-1 α^4
c2 := 0.003 Subscript[M, pl] α m^-1
c3 := 0.0026
ca := 4 (π^2 - 9) * α/(9*\[Pi])


But I get the error

NDSolve::ndsz: At z == 0.00006659917949810859, step size is effectively zero; singularity or stiff system suspected. >>


Following the suggestion of @Michael E2 the integrations runs without errors, however the results is pretty crazy:

LogLogPlot[Y[z] /. sol, {z, 10^-3, 10}, PlotRange -> {{10^-3, 10}, {10^-15, 1}}] • You seem to gave us the wrong code. f[x] is not defined in your code but used in the NDSolve Sep 29 '16 at 14:40
• You are right, I've just edit the post now should be correct. Sep 29 '16 at 14:48
• Okay now a "c3" is undefined for the initial constants. Sep 29 '16 at 14:58
• c3 it's just a constant. Sep 29 '16 at 15:14
• Then give us examples for values for c1, c2, c3, ca please. Sep 29 '16 at 15:46

The relative sizes of the initial values for the derivatives and for the variables, suggests you need (at the least) a high WorkingPrecision, as well as "StiffnessSwitching", taking a hint from the OP's error.

IVs = Rule @@@ odes[[-3 ;;]];
params = {m -> 1000, α -> 137^-1, Subscript[M, pl] -> 2.44*10^18};
{Y'[z], S'[z], T'[z]} /.
First@Solve[odes[[;; 3]], {Y'[z], S'[z], T'[z]}] /.
z -> (5*α^2)/4 /. IVs /. params
(*  {1.76281*10^23, -6.61822*10^26, -1.98517*10^27}  *)

{Y[z], S[z], T[z]} /. z -> (5*α^2)/4 /. IVs /. params
(*  {0.000195865, 1.86648*10^-6, 5.59944*10^-6}  *)


The following integrates without error:

c1 := 2.13100 * Subscript[M, pl] m α^4
c2 := 0.003100 * Subscript[M, pl] α m^5
c3 := 0.0026100;
ca := 4 (π^2 - 9)*α/(9*π)

sol = NDSolve[odes /. params, {Y, S, T}, {z, (5*137^-2)/4, 10},
PrecisionGoal -> 8, AccuracyGoal -> 8, WorkingPrecision -> 32,
Method -> "StiffnessSwitching"];

• Now it runs without errors, but still the results is pretty crazy (see the edit on the question). Sep 30 '16 at 8:37
• @Cervantes Try LogLogPlot[Abs@Y[z] /. sol, {z, 10^-3, 10}, PlotRange -> {{10^-3, 10}, {10^-15, 1}}]: i.stack.imgur.com/MDAFf.png. The logarithm of a negative number is complex and won't be plotted. The solution is oscillatory, which appears as a sequence of humps, in case you think the humps are crazy. Sep 30 '16 at 10:00