I am solving a first-order ODE:$\frac{\partial \rho }{\partial t} = -a \rho^2 + b |A(t)|^2 \rho +c|A(t)|^{2m}$This is the evolution of the plasma density in the presence of a laser pulse (complex amplitude A, which is normalized). The coefficients a, b and c depend on the material and the initial laser intensity. Typically, these coefficients are very large (for example a = 1e13, b=2e16 c=1e7). Now I solve this equation (in Python, and I also played with C) with the conventional Runge-Kutta 4th order method (also tried 6th order); here is a little bit of code (quite standard, nothing fancy):for l in range(Nr):    mag= np.abs(Uo[:,l] )**2     for k in range(1,Nt):          k1 = dT*(b*rho[k-1,l]*mag[k-1] + c*mag[k-1]**m - a*rho[k-1,l **2.0 )           tmp1=0.5  * ( mag[k]+mag[k-1])         tmp2=rho[k-1,l] + 0.5*k1              k2 = dT * (b*tmp2 * tmp1 + c*tmp1**m- a* tmp2**2.0  )                  tmp2=rho[k-1,l] + 0.5*k2            k3 = dT * (b*tmp2 * tmp1 + c*tmp1**m- a* tmp2**2.0 )                   tmp2=rho[k-1,l] + k3                k4 = dT * (b*tmp2 * mag[k] + c*mag[k]**m- a* tmp2**2.0 )           rho[k,l]=   rho[k-1,l] + (k1/6. + k2/3. + k3/3. + k4/6.)I also loop over space (Nr) because I have a time and space dependent plasma density. Since I discretized my time axis (with a total of Nt points, dT is usually 1e-16s), I do not use adaptive step size control (besides, I want this code to be fast).Now, if I do this integration with low power, everything works fine, and I get a nice result (image showing the temporal evolution of plasma density $\rho$ for low power, at one space point):But if I increase the pulse power, my Runge Kutta solution blows up, and I get NaNs. Reason is, that for high powers (say 1e9 Watts), the coefficients become huge (for example a = 1e13, b=2e52,  c=1e13).My question is:Do you have an idea how I can still get a result, without using step size control or using a finer grid? I wouldn't mind using approximations.I am really thankful for any help:)
Runge Kutta solution blows up for a first order ODE with very large coefficients
ode;runge kutta;numerical modelling