Singularity or stiff system in NDSolve

Question

I have the following piece of code

(*UV cuttoff*)
rmax = 10;
(*Fixed input values*)
nc = 4;
nf = 4;
(*Value of α1*)
α1 = 0.95;

b0 = 1/(6 π) (11 nc - 2 (nc + 2) nf) ;

b1 =  1/(24 π^2) (34 nc^2 - 10 nc (nc - 2) nf  -  
     6 (nc + 1) (nc - 2) (nc - 2) nf) ;

eqalp1 = alp1'[mu1] + (b0 alp1[mu1]^2 + b1 alp1[mu1]^3) == 0;


sol1 = NDSolve [{ eqalp1, alp1[0] == α1}, 
   alp1, {mu1, 0, 10^24}];

alph1[mu1_] := First[Evaluate[alp1[mu1] /. sol1]]

alph1[mu1];
α[mu_] := alph1[mu]

It works fine if I choose nf=1,2,3 but when I set nf=4 and higher I get the following error message

I tried to solve the differential equation not from zero but from other points. Every time I get a similar error. I also tried changing the tuned input for α1 in the beginning of the code, but the same message persists.

Can someone suggest something to bypass this?

P.S: I have used different combinations of max steps and working precision.I forgot to mention that previously. Finally, I tried to adopt this answer and use the WhenEvent command but it has no effect whatsoever.

Answer 1

Your ODE can be solved exactly, by the way. The time to reach a singularity is given by the following integral, if it converges:

Integrate[
 -1/(b0 alp1[mu1]^2 + b1 alp1[mu1]^3), 
 {alp1[mu1], α1, Infinity}]

(*  0.451388  *)

Therefore the IVP has a singularity around mu1 = 0.451. This agrees with the numerical integration performed by NDSolve (for nf = 4).

(The integral comes from separating the variables. Given $${d\alpha\over d\mu} = v(\alpha)\,,$$ we get, starting from the initial condition $\alpha(0)=\alpha_1$, $$\int_{\alpha_1}^{\infty} {d\alpha \over v(\alpha)} = \int_0^M d\mu = M$$ The value of $M$, if the integral on the left converges, is the time $\mu = M$ at which $\alpha \rightarrow \infty$.)