Help solving a tricky matrix PDE coupled system I am trying to a solve a matrix equation of the form$$\begin{equation}\label{jacobi_cond}
\frac{\partial \mathbf{P} (x_3)}{\partial x_3} = \mathbf{H}(x_3) \mathbf{P}(x_3).
\end{equation}$$
where
$$\frac{\partial}{\partial x_3}\begin{bmatrix} \chi_1(x_3) & \chi_2(x_3) \\ \chi_3(x_3) & \chi_4(x_3) \end{bmatrix} =\begin{bmatrix} \Delta_1(x_3) & \Delta_2(x_3)  \\ \Delta_3(x_3)  & \Delta_4(x_3) \end{bmatrix}  \begin{bmatrix} \chi_1(x_3) & \chi_2(x_3) \\ \chi_3(x_3) & \chi_4(x_3) \end{bmatrix}$$
with initial condition 
$$\mathbf{P}(x_3=0) = \mathbf{I}$$
the identity matrix. 
The $\chi(x_3)$ entries have real and imaginary parts. The $\Delta(x_3)$ functions are nontrivial. I should note also that $\text{Tr}[\mathbf{H}(x_3)] \neq 0.$
I am trying to get a handle on the determinant of the matrix $\mathbf{P}(x_3).$ 
$\textbf{So far I have tried}$ using Jacobi's formula. This is easy to implement if $\text{Tr}[{\mathbf{H}(x_3)}] = 0,$ meaning the initial condition becomes the determinant, however in my case it does not hold.
I have also tried solving this numerically, but the determinant blows up rapidly. So I am really stuck at this.
I would really appreciate any ideas or discussion. Thanks.
 A: This is actually a linear first-order differential system with variable coefficients, in terms of $x_3 = x$. The present system can be recast as
$
\boldsymbol{\chi}'(x) =  {\bf A}(x)\, \boldsymbol{\chi}(x)
$,
where $\boldsymbol{\chi} = (\chi_1, \chi_2, \chi_3, \chi_4)^\top$ and
$$
{\bf A} = 
\begin{pmatrix}
\Delta_1 & 0 & \Delta_2 & 0\\
0 & \Delta_1 & 0 & \Delta_2\\
\Delta_3 & 0 & \Delta_4 & 0\\
0 & \Delta_3 & 0 & \Delta_4
\end{pmatrix} .
$$
The solutions are given by
$\boldsymbol{\chi}(x) = e^{{\bf B}(x)}\,  \boldsymbol{\chi}_0$,
where $\boldsymbol{\chi}_0 = (1,0,0,1)^\top$ and where $e^{{\bf B}(x)}$ denotes the matrix exponential of ${\bf B}(x) = \int_0^x {\bf A}(\xi)\, \text d \xi$, i.e. the matrix exponential of
$$
{\bf B} = 
\begin{pmatrix}
d_1 & 0 & d_2 & 0\\
0 & d_1 & 0 & d_2\\
d_3 & 0 & d_4 & 0\\
0 & d_3 & 0 & d_4 
\end{pmatrix},
\qquad
d_n(x) = \int_0^x \Delta_n(\xi)\,\text d \xi\, .
$$
Using symbolic calculus software, the following $x$-dependent coefficients $\chi_n$ are obtained:
$$
\boldsymbol{\chi} = \frac{e^{(d_1 + d_4)/2}}{2\sqrt{\mathcal{D}}} \begin{pmatrix}
(d_1-d_4+\sqrt{\mathcal{D}})e^{\sqrt{\mathcal{D}}/2} - (d_1-d_4-\sqrt{\mathcal{D}})e^{-\sqrt{\mathcal{D}}/2} \\
2\, (e^{\sqrt{\mathcal{D}}/2} - e^{-\sqrt{\mathcal{D}}/2})\, d_2 \\
2\, (e^{\sqrt{\mathcal{D}}/2} - e^{-\sqrt{\mathcal{D}}/2})\, d_3 \\
(d_1-d_4+\sqrt{\mathcal{D}})e^{-\sqrt{\mathcal{D}}/2} - (d_1-d_4-\sqrt{\mathcal{D}})e^{\sqrt{\mathcal{D}}/2}
\end{pmatrix}
$$
where $\mathcal{D} = (d_1 - d_4)^2 + 4d_2d_3$. The corresponding Maple code is given below:
with(LinearAlgebra):
B := Matrix(<<d__1,0,d__3,0>|<0,d__1,0,d__3>|<d__2,0,d__4,0>|<0,d__2,0,d__4>>);
M := MatrixExponential(B):
chi := M . Vector(<1,0,0,1>);

