I'm wondering if there's an nice solution to the integral $\int_0^{a} \exp(-s\mathbf{A}) \mathbf{M} \exp(s\mathbf{C}) \: ds$ . The matrices $\mathbf{A}$ and $\mathbf{C}$ are invertible if that matters, but they do not in general commute.
When $\mathbf{A}=\mathbf{C}$, a nice solution was posted here: Derivative of Matrix Exponential as Integral
But I was not able to use to utilize that further.
I had a similar problem lately which I solved by using the "block-triangular method" which came up in a post about the matrix exponential derivative which works for analogous reasons to the calculation of simple matrix exponential integrals discussed in this post.
The solution of your integral can be computed via the upper right block of the matrix exponential:
$$\exp\Bigg(\begin{bmatrix} A & M \\ 0 & B \end{bmatrix} \cdot a \Bigg) = \begin{bmatrix}\exp(a A) & \exp(a A) \int^a_0 \exp(-s A) M \exp(s B) ds \\ 0 & \exp(a B)\end{bmatrix}$$
Why does it work?
The setup of ODEs is slightly different from the original post on the integral of a single matrix exponential:
\begin{align} \dot{X}(t) &= AX(t) + MU(t) \\ \dot{U}(t) &= BU(t) \end{align} $U(t) = \exp(tB)U_0$ is simply the solution of the homogeneous ODE in the augmented variable, which yields a inhomogeneous ODE for $X(t)$:
$$\dot{X}(t) = AX(t) + M\exp(tB)U_0$$,
solved by $X(t) = \exp(t A)X_0 + \exp(t A) \int^t_0 \exp(-\tau A) M \exp(\tau B) d\tau U_0$.
Putting it all back together, the joined variables form a homogeneous ODE:
$$ \begin{bmatrix} \dot{X}(t) \\ \dot{U}(t) \end{bmatrix} = \begin{bmatrix} A & M \\ 0 & B \end{bmatrix}\begin{bmatrix} {X}(t) \\ {U}(t) \end{bmatrix} $$
solved by: $$ \begin{bmatrix} {X}(t) \\ {U}(t) \end{bmatrix} = \exp\Bigg(\begin{bmatrix} A & M \\ 0 & B \end{bmatrix} \cdot t \Bigg) \begin{bmatrix} X_0 \\ U_0 \end{bmatrix} = \begin{bmatrix}\exp(t A) & \exp(t A) \int^t_0 \exp(-\tau A) M \exp(\tau B) d\tau \\ 0 & \exp(t B)\end{bmatrix} \begin{bmatrix} X_0 \\ U_0 \end{bmatrix} $$
In order to compute the integral in question, the final step would be to extract the upper right block of the appropriate matrix exponential and premultiply it with $\exp(-a A)$ to cancel out the unwanted factor in front of the integral.
