I have a system of two differential equations
$$\frac{\partial X(t)}{\partial t}=a_1 A X(t)+a_2X(t) B+a_3 C Y(t)+a_4Y(t) D+a_5$$ $$\frac{\partial Y(t)}{\partial t}=b_1 E X(t)+b_2X(t) F+b_3 G Y(t)+b_4Y(t) H+b_5$$
Where capital letters are quaternions, and small letters are real scalars. Notice that I signify the fact that quaternions are not commutative. How can I solve such a differential equation? I couldn't find literature for that. I'd appreciate someone getting me with these equations on the right foot.
Thanks.
This is an extended version of my comment; it might also approximate being an answer, but probably not.
Both left and right multiplication by a quaternion are $\Bbb{R}$-linear maps on $\Bbb{H}$, and so can be represented by matrices. Let $L_Q$ is the matrix corresponding to left multiplication by some quaternion $Q$, and $R_Q$ the matrix corresponding to right multiplication: that is, for any other quaternion $P$, we should have $L_Q(P)=QP$, $R_Q(P)=PQ$. We can find $L_Q$, $R_Q$ explicitly by noting the effect of multiplication on the basis elements $\{1,i,j,k\}$; if $Q=a+bi+cj+dk$, then $$ Q1=1Q=a+bi+cj+dk $$ and so the first column of both $L_Q,R_Q$ should be $\left(\begin{array}{c}a\\b\\c\\d\end{array}\right)$;
$$ Qi = -b + ai + dj - ck \\ iQ = -b + ai - dj + ck $$ and so the second column of $L_Q$ should be $\left(\begin{array}{c}-b\\a\\-d\\c\end{array}\right)$ while the second column of $R_Q$ is $\left(\begin{array}{c}-b\\a\\d\\-c\end{array}\right)$; and so on.
Proceeding in this way, the eventual result is that $$ L_Q = \left(\begin{array}{cccc} a & -b & -c & -d \\ b & a & -d & c \\ c & d & a & -b \\ d & -c & b & a \end{array}\right), \quad R_Q = \left(\begin{array}{cccc} a & -b & -c & -d \\ b & a & d & -c \\ c & -d & a & b \\ d & c & -b & a \end{array}\right) $$ Note that we are identifying the quaternions being multiplied with column vectors. Since matrices naturally multiply on the left of column vectors, this means $L_Q$ behaves nicely with respect to matrix multiplication but $R_Q$ does not; in fact, we have $L_{Q_1}L_{Q_2}=L_{Q_1Q_2}$, but $R_{Q_1}R_{Q_2}=R_{Q_2Q_1}$. Some sources (e.g., wikipedia) may prefer to identify quaternions with row vectors, in which case all these matrices will be transposed and the opposite things will be nice.
Notice that both $L_Q$ and $R_Q$ are "almost skew-symmetric": that is, they are of the form $aI+S$, where $S=-S^T$. This is an additive property, so in fact any map of the form $L_{Q_1}+R_{Q_2}$ (which corresponds to the quaternionic map $X \mapsto Q_1X+XQ_2$) will also be almost skew-symmetric.
Indeed, any linear map $T$ which is almost skew-symmetric can be written in the form $L_{Q_1}+R_{Q_2}$ for some quaternions $Q_1$, $Q_2$. This isn't too hard to check (basically, the off-diagonal coefficients show up in pairs, and you can pick the sum and difference of each pair), but it's also mostly irrelevant to your question except as an impossibility result (it means that whatever method you come up with to solve your ODE must be general enough to deal with an arbitrary almost skew-symmetric ODE).
I'm not sure how helpful this is for solving your equation, but it certainly would be helpful for solving similar but simpler things. For example, let's say we just had the single quaternionic ODE $X'(t)=AX+XB$. Then you could write this as a vector ODE $X'(t)=(L_A+R_B)X=TX$ for some almost skew-symmetric $T$. This can be solved via the matrix exponential: the general solution for initial condition $X_0$ is $X(t)=e^{Tt}X_0$.
Now, we can write $T=aI+S$ where $S$ is genuinely skew-symmetric; as $aI$ and $S$ commute with each other, we have $e^{Tt}=e^{aIt}e^{St}=e^{at}e^{St}$. It is a well-known fact that the matrix exponential of a skew-symmetric matrix is a rotation matrix (that is, it's orthogonal with determinant $1$), so this means $X(t)$ forms a "logarithmic spiral" in $4$ dimensions: that is, its magnitude varies exponentially with time, and its normalized version is a constant-velocity rotation. (Be careful, however: $4$-dimensional rotations need not be confined to a single plane...)
If we ignore the inhomogeneous terms $a_5,b_5$ in your equation, it corresponds to the block matrix ODE $$ \frac{d}{dt}\left(\begin{array}{c}X\\Y\end{array}\right)=\left(\begin{array}{cc}T_1 & T_2\\T_3 & T_4\end{array}\right)\left(\begin{array}{c}X\\Y\end{array}\right) $$ where the $T_i$ are $4 \times 4$ almost skew-symmetric matrices. I don't know of anything general that can be said about this equation. However, if it turned out that in your specific case $T_2=-T_3$, this would be an $8 \times 8$ almost skew-symmetric matrix, and so its solution would be an $8$-dimensional "logarithmic spiral" of the form described above. This corresponds to the relations $b_1E=-a_3C$, $b_2F=-a_4D$ in your original system.
In either case, after we've solved the homogeneous part of your equation, you can use standard ODE methods — undetermined coefficients or variation of parameters — to re-insert the inhomogeneous terms.