# differential equations – Stiff Sistem of ODEs

I am trying to solve a system of differential equations with NDSolve, but I get as output

NDSolve :: ndsz: in lnP == -21.0018, the step size is effectively zero; singularity or suspicious rigid system.

The core of the program is:

``````DiffEqs := If(Pb0OverP0 != 1, {
lnPb'(lnP) == dlnPbdlnP(lnP),
}, {
(*lnPb'(lnP)(Equal)dlnPbdlnP(lnP),*)
});
CentralValues := If(Pb0OverP0 != 1, {
lnPb(#) == # + Log(Pb0OverP0),
r(#) ==
Sqrt((4.*RhoSung*(10^-8)(*(1-Exp(#))*))/(
3.*(((Epsilon)b(#) + (Epsilon)d(#)) +
PcentralGeo*Exp(#))(*Exp(#)*))),
m(#) ==(*(4. (Pi))/3 *)
r(#)^3/RhoSung *((Epsilon)b(#) + (Epsilon)d(#))
}, {

r(#) ==
Sqrt((3*10^(-8))/(RSg^2*((Epsilon)b(#) +
PcentralGeo*Exp(#))*((Epsilon)b(#) +
3*PcentralGeo*Exp(#)) (-(Lambda) + (4 (Pi))/
PcentralGeo))),
m(#) ==(*(4. (Pi))/3 *)
r(#)^3/RhoSung *((Epsilon)b(#) + (Epsilon)d(#))
}) &;
NDSolveOptions = {MaxSteps -> Infinity,AccuracyGoal -> 10,PrecisionGoal -> 10,
Method -> {"StiffnessSwitching",Method -> {"ExplicitRungeKutta", Automatic}}};
allsol = If(Pb0OverP0 != 1,
NDSolve(Flatten({DiffEqs(#), CentralValues(#)}), {m, r,
lnPb}, {lnP, #, LnPMin}, Sequence@NDSolveOptions),
NDSolve(
Flatten({DiffEqs, CentralValues(#)}), {m, r}, {lnP, #, LnPMin},
Sequence@NDSolveOptions)) & /@ LnP0List;
``````

With

``````Pb0OverP0=0.855924,logPc=-20.4569,PcentralGeo=1.3052*10^-9,(Lambda)=0
LnP0List={0., -0.0921034, -0.184207, -0.27631, -0.368414, -0.460517, -0.55262,
-0.644724, -0.736827, -0.828931, -0.921034, -1.01314, -1.10524,
-1.19734, -1.28945, -1.38155, -1.47365, -1.56576, -1.65786, -1.74996,
-1.84207, -1.93417, -2.02627, -2.11838, -2.21048, -2.30259, -2.39469,
-2.48679, -2.5789, -2.671, -2.7631, -2.85521, -2.94731, -3.03941,
-3.13152, -3.22362, -3.31572, -3.40783, -3.49993, -3.59203, -3.68414,
-3.77624, -3.86834, -3.96045, -4.05255, -4.14465, -4.23676, -4.32886,
-4.42096, -4.51307, -4.60517},LnPMin=-64.918,RhoSung=1.37*10^-8
``````

Essentially, the RhoSung, PcentralGeo and PbarcentralGeo parameters are used to avoid having to deal with too large numbers. logPc is the logarithm of PcentralGeo.
The equations are:

``````    dlnPbdlnP(lnP_) =
Exp(lnP + logPc)/
Exp(lnPb(lnP) + logPbc) (((Epsilon)b(lnP) +
Exp(lnPb(lnP) + logPbc))/(((Epsilon)b(lnP) + (Epsilon)d(lnP)) +
Exp(lnP + logPc)))

Exp(lnP)*3/RhoSung*((Epsilon)b(lnP) + (Epsilon)d(lnP))*
r(lnP)^2/((((Epsilon)b(lnP) + (Epsilon)d(lnP)) +
PcentralGeo*Exp(lnP))/(
2*r(lnP) (m(lnP) - r(lnP)))*(m(lnP)/PcentralGeo + (
3*Exp(lnP)*(r(lnP))^3)/RhoSung) +
2/r(lnP)*(PtFunc(lnP) - Exp(lnP)*PcentralGeo))

lnP)/(3/RhoSung*((Epsilon)b(lnP) + (Epsilon)d(lnP))*r(lnP)^2)
``````

Finally, (epsilon) byd (epsilon) d are defined by interpolating a given table of values. (epsilon) b (lnP) range from 10 ^ -10 to 10 ^ -21 with lnP ranging from 0 to -64, that is

``````(Epsilon)b(lnP_) =
If(Pb0OverP0 != 1, (Epsilon)OfLnPScale(
lnPb(lnP)), (Epsilon)OfLnPScale(lnP))
(Epsilon)d(lnP_) =
If(Pb0OverP0 != 1, (Epsilon)OfLnPScale(
Log(Exp(lnP) - Exp(lnPb(lnP)))), 0)
``````

With other tables, NDSolve works, but for some reason with the previous parameter it does not work.