Solving a linear system of ODEs in matrix form I would like to solve the linear system of ordinary differential equations
$$
\begin{align*}
\dot\mu&=R^\mathrm T\mu \\
\dot\Sigma&=R^\mathrm T\Sigma+\Sigma^\mathrm T R+\mathrm{diag}(R^\mathrm T\mu)
\end{align*}
$$
where $\mu$ is a column vector, $\Sigma$ is a (symmetric) covariance matrix, and $R$ is a square rate matrix. Clearly $\mu(t)=e^{R^\mathrm Tt}\mu(0)$, but is there any way to solve for $\Sigma(t)$?
 A: The matrix equation to solve reads
$$
\dot{\Sigma} = 2\,\mathrm{sym} (R^\mathrm{T} \Sigma) + \Lambda(t) \, ,
$$
where $\Sigma$ is a $n\times n$ real symmetric matrix, $\Lambda(t) = \mathrm{diag}\!\left(R^\mathrm{T}\! \exp\left(t R^\mathrm{T}\right)\mu(0)\right)$ is a diagonal matrix, and $\mathrm{sym}(M) = \frac{1}{2}(M + M^\mathrm{T})$ denotes the symmetric part. One can note that $\dot\Sigma$ is symmetric, independently from the fact that $\Sigma$ is symmetric or not. In the case where $R$ and $\Sigma$ commute, e.g. $R$ is the identity matrix, then the homogeneous solutions are proportional to $\exp\left(2t\,\text{sym}(R)\right)$. In the general case, it seems hopeless to work on the matrix system itself. A strategy consists in rewriting the differential equation in vector form
$$
\dot{V} = A V + b(t) \, ,
$$
where $V$ is a size-$n (n+1)/2$ vector of coordinates of $\Sigma$, $A$ is a square matrix, and $b(t)$ is a vector. To illustrate, let us consider the case of $2\times 2$ matrices
$$
\Sigma = \left(\begin{array}{cc} \sigma_1 & \sigma_2\\ \sigma_2 & \sigma_3\end{array}\right),\quad R = \left(\begin{array}{cc} r_1 & r_2\\ r_3 & r_4\end{array}\right),\quad \Lambda(t) = \left(\begin{array}{cc} \lambda_1(t) & 0\\ 0 & \lambda_2(t)\end{array}\right).
$$
The differential equation in matrix form gives
$$
\frac{d}{dt}
\left(\begin{array}{c}
\sigma_1\\
\sigma_2\\
\sigma_3
\end{array}\right)
=
\left(\begin{array}{ccc}
2 r_1 & 2 r_3 & 0\\
r_2 & r_1 + r_4 & r_3\\
0 & 2 r_2 & 2 r_4
\end{array}\right)\!
\left(\begin{array}{c}
\sigma_1\\
\sigma_2\\
\sigma_3
\end{array}\right)
+
\left(\begin{array}{c}
\lambda_1(t)\\
0\\
\lambda_2(t)
\end{array}\right) ,
$$
which can be solved with classical methods.
A: I have realized that the matrix equation is a continuous differential Lyapunov equation, which is studied in linear-quadratic optimal control. The system above has the following solution:
\begin{align*}
\mu(t)&=e^{R^\text T t}\mu_0 \\
\Sigma(t)&=\int_0^t e^{R^\text T(t-\tau)} \,\text{diag}(R^\text T e^{R^\text T \tau}\mu_0)\,e^{R(t-\tau)}\,d\tau
\end{align*}
assuming the usual initial conditions $\mu(0)=\mu_0$ and $\Sigma(0)=0$.
