Solving the given equation How to find general solution of matrix differential equation $\textbf{X}'(t)=\textbf{A}\textbf{X}(t)+\textbf{B}(t)$ by Green's Function/Matrix Method?
In here, $0<t<t_{max}$
\begin{equation}
\begin{aligned}
&\textbf{X}(t)=\left( \begin{array}{cccc}
x_{1}(t) \\
x_{2}(t) \\
x_{3}(t)\\
x_{4}(t) \end{array} \right), \textbf{A}
=
\left( \begin{array}{cccc}
a_{11} & a_{12} &a_{13} & a_{14}\\
a_{21} & a_{22} & a_{23} & a_{24} \\
a_{31} & a_{32}& a_{33} & a_{34} \\
a_{41} & a_{42} & a_{43} & a_{44}
\end{array} \right),\\
&\textbf{B}(t)=\left( \begin{array}{cccc}
b_{1}(t)\\
b_{2}(t)\\
b_{3}(t)\\
b_{4}(t) \end{array} \right).
\end{aligned}
\end{equation}
initial conditions:
$$x_1(0)=m_1$$
$$x_2(0)=m_2$$
$$x_3(0)=m_3$$
$$x_4(0)=m_4.$$
$\lvert A \rvert >0$
and $A$ has four distinct complex eigenvalues $\lambda_1=c i$,
 $\lambda_2=-c i$,
  $\lambda_3=d i$,  $\lambda_4=-d i$ where $a_{ij},c,d\in \mathbb{R}$ for $i,j=1,2,3,4.$
Respectively, the eigenvectors are $v_k$ for $k=1,2,3,4.$
 A: It is rather simple, however, there will be a Green's matrix (tensor) not a green's function. The solution to the equation 
$$X'-AX=B$$
May be written as
$$X(t)=X_{0}(t)+\int\mathcal{G}(t, t')B(t')dt'$$
Where
$$X'_{0}(t)-AX_{0}(t)=0$$
and $\mathcal{G}(t, t')$ is a $4\times4$ matrix satisfying 
$$\frac{d}{dt}\mathcal{G}(t, t')-A\mathcal{G}(t, t')=I\delta(t-t')$$
Where $I$ is the identity. One way to find $\mathcal{G}(t, t')$ is via the Fourier transform
$$\mathcal{G}(t, t')=\frac{1}{\sqrt{2\pi}}\int_{\mathbb{R}}\hat{\mathcal{G}}(\omega, t')e^{i\omega{t}}d\omega$$
and
$$\delta(t-t')=\frac{1}{2\pi}\int_{\mathbb{R}}e^{i\omega(t-t')}d\omega$$
Which gives
$$[i\omega{I}-A]\hat{\mathcal{G}}(\omega, t')=I\frac{1}{\sqrt{2\pi}}e^{-i\omega{t'}}$$
Hence
$$\hat{\mathcal{G}}(\omega, t')=\frac{1}{\sqrt{2\pi}}[i\omega{I}-A]^{-1}e^{-i\omega{t'}}$$
and
$$\mathcal{G}(t, t')=\frac{1}{2\pi}\int_{\mathbb{R}}[i\omega{I}-A]^{-1}e^{i\omega(t-t')}d\omega$$
Where you also need to invert the matrix $[i\omega{I}-A]$. Another approach would be to diagonalize the matrix $A$, $A=S\Lambda{S}^{-1}$ (which you've already done), then you do the transformation $\mathcal{G}\rightarrow\mathcal{G}'=S^{-1}\mathcal{G}S$ and you obtain
$$\frac{d}{dt}\mathcal{G}'(t, t')-\Lambda{\mathcal{G}'(t, t')}=\delta(t-t')$$
Or component wise 
$$\mathcal{G}'_{lm}(t, t')=X_{0}(t), \ l\neq{m}$$
and
$$\frac{d}{dt}\mathcal{G}'_{ll}(t, t')-\Lambda_{ll}\mathcal{G}'_{ll}(t, t')=\delta(t-t')$$
Which is easy to solve...
