NIntegrate warns about Overflow, Indeterminate, or Infinity

Question

Executing the following code, NDSolve begins to solve a system of two coupled ordinary differential equations for the dependent variables $\rho_0(t)$ and $\lambda_k(t)$. Around $t=-2$ NIntegrate warns that

NIntegrate::inumri: The integrand [...] has evaluated to Overflow, Indeterminate, or Infinity for all sampling points in the region with boundaries {{0.,4.64782*10^14}}.

causing NDSolve to abort. Could some NIntegrate-whiz tell me what's wrong with my integrand? I looked at it in various 3d plots and it appears completely regular. I tried rerunning the calculation with different values for WorkingPrecision and MaxRecursion in case NIntegrate was having difficulties with small numbers, but without success. Any other suggestions how to diagnose this further and get NDSolve to solve my ODEs would be much appreciated!

sqrt[x_?NumericQ,y_?NumericQ]=Piecewise[{{I Sqrt[-x],Re[x]<0&&Re[y]>=0},{-I Sqrt[-x],Re[x]<0&&Re[y]<0}},Sqrt[x]];
Derivative[1,0][sqrt][x_,y_]=1/(2sqrt[x,y]);
Derivative[0,1][sqrt][x_,y_]=0;

αk[1,1][m2_,g2_,z_,c_]:=1/2 (1/c+m2/z-I g2/z)+I sqrt[1/(c z)-1/4 (1/c-m2/z+I g2/z)^2,-(1/c-m2/z)];
αk[1,2][m2_,g2_,z_,c_]:=1/2 (1/c+m2/z+I g2/z)-I sqrt[1/(c z)-1/4 (1/c-m2/z-I g2/z)^2,(1/c-m2/z)];
αk[2,1][m2_,g2_,z_,c_]:=1/2 (1/c+m2/z-I g2/z)-I sqrt[1/(c z)-1/4 (1/c-m2/z+I g2/z)^2,-(1/c-m2/z)];
αk[2,2][m2_,g2_,z_,c_]:=1/2 (1/c+m2/z+I g2/z)+I sqrt[1/(c z)-1/4 (1/c-m2/z-I g2/z)^2,(1/c-m2/z)];

w={1,1,1,1,-1,-1,-1,-1,-2,2};
μ=Join[Flatten@Table[αk[i,j][m2/s^2,g2/s^2,z,c]s^2,{s,{k,Λ}},{i,2},{j,2}],{k^2/c,Λ^2/c}];

ΔkerT=T Sum[w[[j]]Log[Exp[1/T Sqrt[p^2+μ[[j]]]]-1],{j,10}];
kerT[0][m2_,g2_,z_,c_,T_,p_]=k D[ΔkerT,k]/.k->1;
kerT[j_][m2_,g2_,z_,c_,T_,p_]:=(-1)^j Derivative[j,0,0,0,0,0][kerT[0]][m2,g2,z,c,T,p];

tfI[j_][m2_,g2_,z_,T_,c_,d_,pmax_,ops___]:=(2π^(d/2))/((2π)^d Gamma[d/2]) NIntegrate[kerT[j][m2,g2,z,c,T,p] p^(d-1),{p,0,pmax},ops];

ρ0Flow[ρ0_?NumericQ,λk_?NumericQ,ηk_?NumericQ,g2_?NumericQ,Z1_?NumericQ,T_?NumericQ,c_?NumericQ,d_?NumericQ,N_?NumericQ,pmax_,ops___]:=-(2+ηk)ρ0+(3/2 tfI[1][2ρ0 λk,g2,Z1,T,c,d,pmax,ops]+(N-1)/2tfI[1][0,0,1,T,c,d,pmax,ops])
λkFlow[ρ0_?NumericQ,λk_?NumericQ,ηk_?NumericQ,g2_?NumericQ,Z1_?NumericQ,T_?NumericQ,c_?NumericQ,d_?NumericQ,N_?NumericQ,pmax_,ops___]:=2ηk λk+λk^2 (9/2 tfI[2][2ρ0 λk,g2,Z1,T,c,d,pmax,ops]+(N-1)/2tfI[2][0,0,1,T,c,d,pmax,ops])

AbsoluteTiming@Module[{runner = 0, counter = 0},
  With[{ηk = 0, g2 = 0, Z1 = 1, T = 0.1 E^t, c = 1, d = 3, N = 2,
     pmax = ∞}, 
   run[1] = 
    NDSolve[{ρ0'[
        t] == ρ0Flow[ρ0[t], λk[t], ηk, g2, Z1, 
        T, c, d, N, pmax], λk'[
        t] == λkFlow[ρ0[t], λk[t], ηk, g2, 
        Z1, T, c, d, N, pmax], ρ0[0] == 0.02, λk[0] == 
       0.5}, {ρ0, λk}, {t, -10, 0}, 
     StepMonitor :> 
      counter++ If[Abs[t] > runner, 
        Print@Chop[{counter, Round[t, 1], 
           E^(2 t) ρ0[t], λk[t]}]; runner++]]]]

Answer 1

NIntegrate fails when the arguments of exponentials in the integrand become too large, causing overflows. This can be seen by removing unnecessary code from AbsoluteTiming@NDSolve ... and enclosing it in EvaluationData as follows.

With[{ηk = 0, g2 = 0, Z1 = 1, T = 0.1 E^t, c = 1, d = 3, n = 2, pmax = ∞},  
    EvaluationData[NDSolve[{ρ0'[t] == ρ0Flow[ρ0[t], λk[t], ηk, g2, Z1, T, c, d, n, pmax], 
    λk'[t] == λkFlow[ρ0[t], λk[t], ηk, g2, Z1, T, c, d, n, pmax], 
    ρ0[0] == 0.02, λk[0] == 0.5}, {ρ0, λk}, {t, -10, 0}]]];

Which produces an enormous expression and many error messages. (Note that N has been replaced by n everywhere in the code in the question, because N is a Mathematica function. Using this reserved symbol is not, however, the source of the problem in the question.) The NIntegrate integrand that caused the error then is obtained from

Cases[%, Hold[Message[NIntegrate::inumri, z__]] -> z, Infinity, 1]

which also is very large but contains 18 terms of the general form

-(((29.0983 + 50.3998 I) E^(116.393 Sqrt[1/2 - (I Sqrt[3])/2 + p^2]))/
    ((-1 + E^(116.393 Sqrt[1/2 - (I Sqrt[3])/2 + p^2])) (1/2 - (I Sqrt[3])/2 + p^2)^(5/2)))

Moreover, the limits of integration over p are {0., 4.64782*10^14}. Thus, the exponentials overflow. Probably, the solution is to restructure the integrand so that it consists instead of terms of the general form

-(((29.0983 + 50.3998 I)/
   ((-E^(-116.393 Sqrt[1/2 - (I Sqrt[3])/2 + p^2]) + 1) (1/2 - (I Sqrt[3])/2 + p^2)^(5/2)))

which should avoid the overflows. I hope that this forensic analysis it helpful and regret that I do not have time to pursue the problem futher.