Background: A discrete-time equation that arises in control is apparently the matrix equation: $$ U + A^TUA = - C, $$ where $A$ and $C$ are given real square matrices and the variable is $U$. It is known that if the spectral radius $\rho(A) < 1$, then there is a unique solution $U = \sum_{j \geq 0} (A^j)^T C (A^j)^T$.
Seeing how this arises can be done in many ways. For example, one looks at the case when the dimension is $1$, $u + a^2 u = -c$, and so $u = -(1 -a^2)^{-1} c = -\sum_{j \geq 0} a^jca^j$, using the formula for a geometric series and that $\rho(a) = |a| < 1$, and so then one can analogize and verify that the ansatz $\sum_{j \geq 0} (A^j)^T C (A^j)^T$. Yet another way to see how this arises is to conjugate the matrix equation itself by $A$, which says: $$ A^T U A + (A^2)^TUA^2 = -A^T C A, $$ and by adding the two equations, and inducting, we arrive at $U = \sum_{j=0}^N (A^j)^T C (A^j)^T + (A^{N+1})^TU A^{N+1}$. The second term necessarily tends to 0 with $N \to \infty$, since $\rho(A) < 1$ and by Gelfand's formula. And so then we see that the given formula arises naturally. (The uniqueness can be done by rewriting the equation as a linear system in the $n^2$ variables $\mathrm{vec}(U)$.
Question: I would like to understand how one can make similar manipulations to determine the solution of the continuous-time Lyapunov equation, $$ AU + UA^T = -C, $$ where again $A$ and $C$ are given real square matrices, $C$ is additionally symmetric, and the variable is $U$. It is apparently the case that $$ U = \int_0^\infty e^{At} C e^{A^T t} \, dt, $$ is a solution, provided that the eigenvalues of $A$ lie in the left open complex plane (also called Hurwitz). But is there an easy way to see this is true? Specifically, I would also like to have a more intuitive understanding this solution, for example, ideally I would like to understand why it is natural to guess that this choice of $U$ will satisfy the equation (as opposed to assuming this is the solution and then verifying it).
Start from the Lyapunov equation $$AU+UA^T = -C.$$ Next, left multiply by $e^{At}$ and right multiply by $e^{A^Tt}$ to get $$e^{At}AUe^{A^Tt}+e^{At}UA^Te^{A^Tt} = -e^{At}Ce^{A^Tt}.$$ Since $\dfrac{d}{dt}e^{At} = e^{At}A$ and $\dfrac{d}{dt}e^{A^Tt} = A^Te^{A^Tt}$, this can be written as $$\dfrac{d}{dt}\left[e^{At}\right]Ue^{A^Tt}+e^{At}U\dfrac{d}{dt}\left[e^{A^Tt}\right] = -e^{At}Ce^{A^Tt}.$$ Using the product rule, this can be rewritten as $$\dfrac{d}{dt}\left[e^{At}Ue^{A^Tt}\right] = -e^{At}Ce^{A^Tt}.$$ Finally, integrate from $t = 0$ to $\infty$ to get $$\int_{0}^{\infty}\dfrac{d}{dt}\left[e^{At}Ue^{A^Tt}\right]\,dt = -\int_{0}^{\infty}e^{At}Ce^{A^Tt}\,dt,$$ which simplifies to $$U = \int_{0}^{\infty}e^{At}Ce^{A^Tt}\,dt.$$ Note that this assumes that $e^{At} \to 0$ and $e^{A^Tt} \to 0$ as $t \to \infty$, which only holds when all the eigenvalues of $A$ lie in the left open complex plane.