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),
m'(lnP) == dmadimdlnP(lnP),
r'(lnP) == dradimdlnP(lnP)
}, {
(*lnPb'(lnP)(Equal)dlnPbdlnP(lnP),*)
m'(lnP) == dmadimdlnP(lnP),
r'(lnP) == dradimdlnP(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)))  

dmadimdlnP(lnP_) = 
 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))

dradimdlnP(lnP_) = 
 dmadimdlnP(
   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.