Solution of a first- and second-order consecutive reaction scheme I have a compartmental model system X, a simplified version of which, looks like this:

$$
\text{System }x: \enspace \text{STO} \xrightleftharpoons[k_2]{k_1}  \text{blood} \xrightarrow[]{k_3} \text{urine} 
\equiv
x_1 \xrightleftharpoons[k_2]{k_1} x_2 \xrightarrow[]{k_3} x_3
$$

The system of differential equations is then:
$$
\frac{dx_1}{dt} = -k_1 x_1 + k_2 x_2 \\
\frac{dx_2}{dt} = k_1 x_1 - k_2 x_2 - k_3 x_2\\
\frac{dx_3}{dt} = k_3 x_2 \\
$$

I can solve this system of equations using the following matrix approach.

First, I write the rate matrix [R]. From [R] one can obtain a new matrix [A] by first replacing each diagonal element of [R] by the negative of the sum of each of row elements, and then transposing it:
$$
[R] = 
\begin{bmatrix}
    0            &amp; k_2    &amp; 0 \\
    k_1          &amp; 0      &amp; k_3 \\
    0            &amp; 0      &amp; 0 
\end{bmatrix} \\
[A'] = 
\begin{bmatrix}
    -k_2            &amp; k_2    &amp; 0 \\
    k_1          &amp; -(k_1 + k_3)      &amp; k_3 \\
    0            &amp; 0      &amp; 0 
\end{bmatrix} \\
[A] = 
\begin{bmatrix}
    -k_2            &amp; k_1    &amp; 0 \\
    k_1          &amp; -(k_1 + k_3)      &amp; 0 \\
    0            &amp; k_3      &amp; 0 
\end{bmatrix} \\
$$

I can calculate the amount in each compartment by doing the following:

$$
x(t) = e^{[A]t}x(0) \enspace \text{where} \enspace x(0) = \begin{bmatrix} 0 \\ 1 \\0 
 \end{bmatrix} 
$$

In python:

RMatrix = model_matrix.as_matrix()
row, col = np.diag_indices_from(RMatrix)
RMatrix[row, col] = -(RMatrix.sum(axis=1)-RMatrix[row,col])
AMatrix = RMatrix.T

def content(t):
    cont = np.dot(linalg.expm(t*AMatrix), x0))


This method is working well for me.

Now, I have a little more complicated model where reactants in compartments 1 and 2 of Systems X and Y combine to get product in System Z. 

$X + Y \rightarrow Z$, with a reaction constant of $k_R$.

$$
\begin{align}
&amp;x_1 &amp; \xrightleftharpoons[k_2]{k_1} \enspace &amp; x_2 &amp;\xrightarrow[]{k_3} &amp; x_3 \\
&amp;+  &amp; \ &amp;  + \\
&amp;y_1 &amp; \xrightleftharpoons[k_5]{k_4} \enspace &amp; y_2 &amp;\xrightarrow[]{k_6} &amp; y_3 \\
&amp;\downarrow^{k_R} &amp; \  &amp; \downarrow^{k_R} \\
&amp;z_1 &amp; \xrightleftharpoons[k_8]{k_7} \enspace  &amp; z_2 &amp;\xrightarrow[]{k_9} &amp; z_3 \\
\end{align}
$$

, and the corresponding system of differential equations would be:

$$
\begin{align}
\frac{dx_1}{dt} &amp;= -k_1 x_1 + k_2 x_2 - k_R x_1 y_1 \\
\frac{dx_2}{dt} &amp;= k_1 x_1 - k_2 x_2 - k_3 x_2 -  k_R x_2 y_2\\
\frac{dx_3}{dt} &amp;= k_3 x_2 \\ 
\frac{dy_1}{dt} &amp;= -k_4 y_1 + k_5 y_2 - k_R x_1 y_1 \\
\vdots&amp; \\ 
\frac{dz_z}{dt} &amp;= -k_7 z_1 + k_8 z_2 + k_R x_1 y_1 \\
\end{align}
$$

I am struggling with a method to solve this system of differential equations (1st and 2nd order) to calculate the amount in each compartment at a certain time t, given the initial conditions, $k_R$, and the transfer rates $k_1$, $k_2$, $k_3$, etc...

Can I solve it using the matrix method like the one above for a system of first order differential equations? What other options in Python do I have?

Thanks in advance!
Look for equation solvers in Numpy/ Scipy and load at the top of your script scipy.integrate import odeint and import numpy as np . There are examples of how to numerically solve differential equations in the examples and on line. The Master equation approach does not work for second order steps. There is also the Gillespie method which is a sort of Monte-Carlo method see this answer which has an example code Probability of the single molecule reaction to happen