We have the following cascade of unimolecular chemical reactions
$$\ce{X_1} \xrightarrow{k_1} \ce{X}_2 \xrightarrow{k_2} \ce{X}_3 \xrightarrow{k_3} \ce{X}_4$$
where $k_1, k_2, k_3 > 0$ are the (distinct) rate constants. Assuming that every chemical reaction in this cascade has mass action kinetics, then we have the following system of linear ODEs
$$\begin{bmatrix} \dot x_1\ \dot x_2\ \dot x_3\ \dot x_4\end{bmatrix} = \begin{bmatrix} - k_1 & 0 & 0 & 0\ k_1 & - k_2 & 0 & 0\ 0 & k_2 & - k_3 & 0\ 0 & 0 & k_3 & 0\end{bmatrix} \begin{bmatrix} x_1\ x_2\ x_3\ x_4\end{bmatrix}$$
where $x_i := [\ce{X}i]$ is the (time-varying) concentration of the $i$-th species. Let the initial concentrations be $x_1 (0) =: x{10} \neq 0$ and $x_2 (0) = x_3 (0) = x_4 (0) = 0$. Integrating the linear ODEs, we obtain
$$\begin{array}{rl} x_1 (t) &= x_{10} \, e^{- k_{1} t}\ x_2 (t) &= x_{10} \left(\frac{k_{1}}{k_{2} - k_{1}}\right) e^{- k_{1} t} + x_{10} \left(\frac{k_{1}}{k_{1} - k_{2}}\right) e^{- k_{2} t}\ x_3 (t) &= x_{10} \left( \frac{k_{1} k_{2}}{ \left(k_{2} - k_{1}\right) \left(k_{3} - k_{1}\right)}\right) e^{- k_{1} t} + x_{10} \left( \frac{k_{1} k_{2}}{ \left(k_{1} - k_{2}\right) \left(k_{3} - k_{2}\right)}\right) e^{- k_{2} t} + x_{10} \left( \frac{k_{1} k_{2}}{ \left(k_{1} - k_{3}\right) \left(k_{2} - k_{3}\right)}\right) e^{- k_{3} t}\ x_4 (t) &= x_{10} \left( \frac{k_{2} k_{3}}{ \left(k_{2} - k_{1}\right) \left(k_{3} - k_{1}\right)}\right) \left( 1 - e^{- k_{1} t} \right) + x_{10} \left( \frac{k_{1} k_{3}}{ \left(k_{1} - k_{2}\right) \left(k_{3} - k_{2}\right)}\right) \left( 1 - e^{- k_{2} t} \right) + x_{10} \left( \frac{k_{1} k_{2}}{ \left(k_{1} - k_{3}\right) \left(k_{2} - k_{3}\right)}\right) \left( 1 - e^{- k_{3} t} \right)\end{array}$$
Python code
The following SymPy script
from sympy import *
t = Symbol('t')
x10 = Symbol('x10')
k1, k2, k3 = symbols('k1 k2 k3')
A = Matrix([[-k1, 0, 0, 0],
[ k1,-k2, 0, 0],
[ 0, k2,-k3, 0],
[ 0, 0, k3, 0]])
V, D = A. 374