Solving rate equations with RK method. Don't understand how.

[ale_yoman]

I can't understand how authors of the article solve the system of rate equations with Runge-Kutta method.

Pulsed laser light is transmited through nonlinear absorbing sample. So, here is the system of rate equations which governs the populations of the levels in the three-level system.

$\frac{dN_{0}}{dt} = -\frac{\sigma_{0} I N_{0}}{h \omega} - \frac {\beta I^{2} } {2 h \omega} - \frac{\gamma I^{3} }{3 h \omega} + \frac{N_{1}}{\tau_{s1}}$

$\frac{dN_{1}}{dt} = \frac{\sigma_{0} I N_{0}}{h \omega} - \frac{\sigma_{1} I N_{1}}{h \omega} - \frac{N_{1}}{\tau_{s1}} +\frac{N_{2}}{\tau_{sn}}$

$\frac{dN_{2}}{dt} = \frac{\sigma_{1} I N_{1}}{h \omega} + \frac {\beta I^{2} } {2 h \omega} + \frac{\gamma I^{3} }{3 h \omega} - \frac{N_{2}}{\tau_{sn}}$

Intensity transmitted through the sample is given by
$\frac{dI}{dz} = -\sigma_{0} I N_{0} - \sigma_{1} I N_{1} - \beta I^{2} - \gamma I^{3}$

where $I(r,z,t) = I_{00} \left( \frac{w_{0}^{2}}{w^{2}(z)} \right) \ \exp \left( - \frac{t^{2}}{\tau_{p}^{2}} \right) \ \exp \left( - \frac{2r^{2}}{w^{2}(z)} \right)$

Here $\sigma_{i}$ is the cross-section of the corresponding state, $\beta$ and $\gamma$ - two- and three-photon absoption coefficients, $w(z)$ - radius of the laser (assumed to be Gaussian) beam, $N_{i}(t)$ - population of the i-th state, $\tau_{i}$ - life-time of the corresponding state.

Here we come to the point where my incomprehension is rising. It is stated in the article that the differential equations were first decoupled and then integrated over time, length, and along the radial direction before solving them numerically with Runge-Kutta fourth order method. Limits of integration of $r, t, z$ are $(0,\infty), (-\infty, \infty)$ and $(0,L)$,respectively, where L is the sample length.

So, i wonder how i should do this "de-coupling"?

I thought that the RK method is aplicable for the systems where the derivatives are taken with respect to the same variable. But in this exaple the state population is the function of time, and in the equation for the intensity there is the coordinate derivative. May be I'm wrong?

 

[I like Serena]

Hi ale_yoman! Welcome to PF!

To apply RK you need to uncouple your equations, so they become of the form:
$$y' = f(t, y)$$

Or more specifically:
$$\frac d {dt} \begin{pmatrix}y_1 \\ y_2 \\ ... \\ y_n \end{pmatrix} = \vec f(t, y_1, y_2, ..., y_n)$$

Can you see how this matches your set of equations?

Next the actual application of RK is to start with an initial position y, and step t up with a step size, find the corresponding path that y takes.