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)$?

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.
assuming the usual initial conditions $\mu(0)=\mu_0$ and $\Sigma(0)=0$.