I want to numerically integrate the following differential equation of a covariance matrix $C$,
$$ \dot{C}(t) = J(t)^TC(t) + CJ(t) + \Sigma(t), $$ with initial condition $C(0) = 0$. Here $\Sigma(t)$ is a covariance matrix as well, and $J$ is a Jacobian matrix. The matrices can be assumed to have a maximum dimension of 5$\times$5 for now, and even $\Sigma(t) = \Sigma$.
Are there any standard numerical methods for this problem? Since I'm only interested in the full covariance at the final time $T$, i.e. I only need to know $C(T)$, I thought about integrating only the upper triangular part of $C(t)$, and then at time $T$, multiply the upper triangular part and it's transpose to obtain $C(T)$. This will make the matrix symmetric, however, I'm not sure that it will preserve positive semi-definiteness.