NDSolve Step Size Effectively Zero

Question

I'm trying to solve a PDE for a function of 2 variables. The most accurate parameterizations of this equation are very unwieldy and involve numerous piecewise elements, and so right now I'm trying to solve an approximate form. I'm getting step size effectively 0 errors, unless I increase the starting step size, in which case the solution I get is clearly unphysical. Since this is an astrophysics problem, I'm working over many orders of magnitude, which I suspect might be part of the issue.

Code:

z[x_]:= 2458.31 - 100.087 x + 1.23213 x^2 - 0.0046743 x^3

lamavg[t_] := Min[1, 0.01 + 0.07 z[t]];

Ufit2[M_,t_]:=0.50519 + 3.127*10^10/M^2 - 274337./M + 2.12127*10^-10 M           
-1.92858*10^-20 M^2 - 6.20762*10^-11 t;

fixedM[x_] := 
0.00006192808740866853/
591.5967816994163 (32.59434080693661 - 2.0297457454952188*^-7 x + 
1.2887918124478197*^-15 x^2 - 2.5861699533219344*^-25 x^3 + 
1.941139394441828*^-35 x^4 - 5.134757798851362*^-46 x^5)

DifEq = D[P[M, t],t] == -M l /sol^2 D[(1 - e)/e                                  
lamavg[t] *3.93242 * Ufit2[M, t] P[M, t], M] //. {e -> 
0.1,l -> 1.26 * 10^(31), sol -> 3*^8 };

soln = NDSolve[{DifEq, P[M, 8.15240949872944*^8] == 
1.629652976233948/M, P[1*^5 .263, t] == fixedM[t]}, P[M,
t], {t, 1000, 8.15240949872944*^8}, {M, 1*10^5.263, 10000000000}, 
Method -> { "StiffnessSwitching", "NonstiffTest" -> False}]

In terms of what I expect, this is a black hole mass function over the evolution of the universe, so P shouldn't fall off as time decreases, and there should be a difference of several orders of magnitude (especially at large t) between the low and high mass end.

Also I shortened the interval because the part here is the most important part, but ideally I would be able to continue the integration up to t of about 10 billion (not important at the moment though).

Thanks!

Answer 1

This evaluates without errors:

z[x_] := 2458.31 - 100.087 x + 1.23213 x^2 - 0.0046743 x^3

lamavg[t_] := Min[1, 0.01 + 0.07 z[t]]

Ufit2[M_, t_] := 0.50519 + 3.127*10^10/M^2 - 274337/M + 2.12127*10^-10 M - 
                 1.92858*10^-20 M^2 - 6.20762*10^-11 t;

e = 1/100;
l = (126/100)*10^(31);
sol = 3*^8;

DifEq = D[P[M, t],t] == -M l/ sol^2 D[(1 - e)/e lamavg[t]*3.93242*Ufit2[M, t] P[M, t], M];

fixedM[x_] := 0.00006192808740866853/ 591.5967816994163 
   (32.59434080693661 - 2.0297457454952188*^-7 x + 
    1.2887918124478197*^-15 x^2 - 2.5861699533219344*^-25 x^3 + 
    1.941139394441828*^-35 x^4 - 5.134757798851362*^-46 x^5)

soln = NDSolve[{
   DifEq,
   P[M, 8.15240949872944*^8] == 1.629652976233948/M,
   P[0.263*^5, t] == fixedM[t]},
  P[M, t], {M, 0.263*^5, 10000000000}, {t, 1000, 8.15240949872944*^8},
   Method -> "StiffnessSwitching"]

Plot3D[P[M, t] /. soln, {M, 0.263*^5, 10000000000}, {t, 1000, 8.15240949872944*^8}]