I see that Matlab uses the spectral decomposition to solve the continuous Lyapunov equation
$$AX+XA=B$$
The formula they use for positive definite $A$ with matrix of eigenvectors $U$ and column vector of eigenvalues $s$ is as follows
$$U \left( \frac{U' BU}{s + s'} \right) U'$$
(Addition of row and column is done with NumPy broadcasting rules)
Any suggestions on how to derive this?