A set of three differential equations is described as
\begin{equation} \frac{\mathrm d}{\mathrm dt}\ \begin{pmatrix} X(t)\\ Y(t)\\ Z(t) \end{pmatrix} = \begin{pmatrix} -a-bS & -ibD &0\\ -ibD & -a+bS & C\\ 0 & -C & -2a \end{pmatrix} \begin{pmatrix} X(t)\\Y(t)\\Z(t) \end{pmatrix} + \begin{pmatrix} 0 \\ 0\\ -2b \end{pmatrix} \end{equation} Here, $a,b,C,D$ and $S$ are independent of $t$. How can one approach to solve these equations with $general$ initial conditons $X(0)=X_0$, $Y(0)=Y_0$ and $Z(0)=Z_0$?
When you have a linear ODE-system
$$\dot{X} = AX + b,$$
(with $A \in \mathbb{C}^{n \times n},\ X,b \in \mathbb{C}^n$) it is a common strategy to look first at the associated homogenous system
$$\dot{X} = AX.$$
You can make an ansatz for the general solution in form of
$$X_{\text{hom}} = \mathrm e^{tA} C,$$
whereby $C \in \mathbb{C}^n$. For computing $\mathrm e^{tA}$, transform $A$ in Jordan canonical form
$$A = P^{-1} J P,$$
and note that every $J$ can be uniquely written as a sum of a diagonal matrix $D$ and a nilpotent matrix $N$:
$$ J = D + N.$$
Since $D$ and $N$ commute, one can easily show
$$\mathrm e^{tA} = P^{-1} \mathrm{e}^{tD} \mathrm{e}^{tN} P,$$
which can be calculated.
Now, for the inhomogeneous system one shows with variation of constants
$$\dot{C} = (\mathrm e^{tA})^{-1} b,$$
so
$$ C = \int (\mathrm e^{tA})^{-1} b + K) \ \mathrm d t,$$
whereby $K \in \mathbb{C}^{n \times n}$. Now
$$X = \mathrm e^{tA} \left (\int (\mathrm e^{tA})^{-1} b + K) \ \mathrm d t\right).$$
You now can plug in the initial values in order to get $K$.