I have the following coupled linear DE that I wish to solve. However, their coefficients are non-constant, namely, I wish to solve
$$ \frac{d}{dx}\begin{bmatrix}y(x) \newline z(x)\end{bmatrix} = \begin{bmatrix}C_{1}+C_{2}x^{2} & C_{3} \newline D_{3} & D_{1}+D_{2}x^{2}\end{bmatrix}\begin{bmatrix}y(x) \newline z(x)\end{bmatrix} $$
where $C_{j}$ and $D_{j}$ are constants (not $x$ dependent). I tried solving it the old fashion way by finding the eigenvalues and eigenvectors of the matrix, but they became complicated and seems intractable. I was thinking about Laplace transforms but I would have to do two Laplace transforms for the $y(x)$ and $z(x)$ equations respectively and that doesn't seem like it will be useful. How should I go about solving the preceding equations?
I don't expect to find any 'nice' solutions for general values of the constants. Note that in the particular case $C_1=C_2=0$ and $C_3=1$, the system is equivalent to the second order ODE $$y''-(D_1+D_2x^2)y'-D_3y=0$$ which is a case of an Triconfluent Heun equation. Note that it's solutions (called Heun functions) are, in some sense, more complicated than Hypergeometric functions as their group of symmetries is of order 196 (in opposite of 24) and their coefficients in the power series solutions do not obeys simple two-term recurrence relations (in opposite of what happens with the larger class of Generalized hypergeometric functions).