In short, I'm trying to find some method, apart from numerically integrating, to find the value of \begin{equation} \bar{X}_t = \int_0^t e^{A^T (t-s)} Q e^{As} ds, \end{equation} where $Q$ is symmetric and positive definite, and where $A$ is a stable matrix. (That is, eigenvalues all below zero.)
Background knowledge
I already know that $X$, which I defined as \begin{equation} X = \int_0^\infty e^{A^T (t-s)} Q e^{A (t-s)} ds = \int_0^\infty e^{A^T s} Q e^{A s} ds, \end{equation} can be found by solving the Lyapunov equation \begin{equation} A^T X + X A + Q = 0. \end{equation} Similarly, I also know that $X_t$, defined as \begin{equation} X_t = \int_0^t e^{A^T s} Q e^{A s} ds, \end{equation} can be found through \begin{equation} X_t = X - e^{A^T t} X e^{A t} \end{equation} or alternatively through \begin{equation} A^T X_t + X_t A + Q - e^{A^T t} Q e^{A t} = 0. \end{equation} But after several days of searching, I haven't found a method of finding $\bar{X}_t$. Now I'm curious: is there a method (apart from numerical integration) to find $\bar{X}_t$?
It took a few more days, but I figured it out.
First, we want to find a matrix of eigenvectors $V$ and corresponding eigenvalues $D$ such that $A = VDV^{-1}$. (This can be done with the Matlab 'eig' function.) Then we have \begin{equation} e^{At} = e^{VDV^{-1}t} = Ve^{Dt}V^{-1}. \end{equation} This turns the integral into \begin{eqnarray} \bar{X}_t & = & \int_0^t e^{A^T (t-s)} Q e^{As} \, ds \\ & = & \int_0^t e^{(VDV^{-1})^T (t-s)} Q e^{VDV^{-1} s} \, ds \\ & = & V^{-T} \left(\int_0^t e^{D^T (t-s)} V^T Q V e^{D s} \, ds\right) V^{-1}. \end{eqnarray} Note that $D$ is a symmetric matrix. However, we are using the conjugate transpose here, so we do not have $D = D^T$. After all, eigenvalues may be complex.
The above integral may be solved element-wise. If we define \begin{eqnarray} \bar{Q} & = & V^T Q V, \\ \tilde{Q} & = & \int_0^t e^{D^T (t-s)} V^T Q V e^{D s} \, ds, \end{eqnarray} then the elements of $\tilde{Q}$ can be found from the elements of $\bar{Q}$ through \begin{equation} \tilde{Q}_{i,j} = \begin{cases} \bar{Q}_{i,j} \left(\frac{e^{\bar{\lambda}_i t} - e^{\lambda_j t}}{\bar{\lambda}_i - \lambda_j}\right) & \mbox{if } \bar{\lambda}_i \neq \lambda_j, \\ \bar{Q}_{i,j} e^{\lambda_j T} T & \mbox{if } \bar{\lambda}_i = \lambda_j. \end{cases} \end{equation} Here, the $\lambda$ values are the eigenvalues of $A$, and hence form the diagonal of $D$. Also, $\bar{\lambda}$ denotes the complex conjugate of $\lambda$.
Once $\tilde{Q}$ has been found, the final solution is \begin{equation} \bar{X}_t = V^{-T} \tilde{Q} V^{-1}. \end{equation}