How to solve an ODE of this form $\dot{P}(t)=A^TP(t)+P(t)A$? Background
If $\dot{x}(t)=A \,x(t)$, then we know the solution is $x(t)=e^{At}x(0)$. 
Question
Now let $\dot{P}(t)=A^TP(t)+P(t)A$, what is $P(t)$?
Attempt
If we take $P(t)$ as common factor (I'm not sure I'm doing this correctly though) ,then we have: 
$$
\dot{P}(t)=A^TP(t)+P(t)A=(A^T+A)P(t)
$$ 
and using the same rationale in the background section, we have 
$$
P(t)=e^{(A^T+A) t}P(0)=e^{A^T t}P(0)e^{At}
$$
Is this a correct solution?
Note
I know the answer is $P(t)=e^{A^T t}P(0)e^{At}$ but I'm not sure how to find it.
 A: Consider the differential equation
$$
\dot{P}(t)=A^TP(t)+P(t)A
$$
We can rewrite it as:
$$
\dot{P}(\tau)-A^TP(\tau)-P(\tau)A=0
$$
Multiply the last differential equation by $e^{-A^T\tau}$ from the left and by $e^{-A\tau}$ from the right then 
\begin{equation}
e^{-A^T\tau}\dot{P}(\tau)e^{-A\tau}-e^{-A^T\tau}A^TP(\tau)e^{-A\tau}-e^{-A^T\tau}P(\tau)Ae^{-A\tau}=0 \quad (1)
\end{equation}
We can easily see that equation $(1)$ is generated by the derivative of three multiplied functions:
 \begin{equation}
\frac{d}{d\tau}\bigg(e^{-A^T\tau}{P}(\tau)e^{-A\tau}\bigg)=e^{-A^T\tau}\dot{P}(\tau)e^{-A\tau}-A^Te^{-A^T\tau}P(\tau)e^{-A\tau}-e^{-A^T\tau}P(\tau)Ae^{-A\tau}=0
\end{equation}
So we come up with:
$0=\frac{d}{d\tau}\bigg(e^{-A^T\tau}{P}(\tau)e^{-A\tau}\bigg)$
Take the integral of both sides
\begin{equation}
0=\int_0^t\frac{d}{d\tau}\bigg(e^{-A^T\tau}{P}(\tau)e^{-A\tau}\bigg)d\tau=\bigg(e^{-A^T\tau}{P}(\tau)e^{-A\tau}\bigg)\bigg]_0^t=e^{-A^Tt}{P} \tau)e^{-At}-e^{0}{P}(0)e^{0}\end{equation}
 \
\begin{equation}
0=e^{-A^Tt}{P}(t)e^{-At}-P(0)
\end{equation}
Multiply left and right sides by $e^{A^Tt}$ and $e^{At}$ respectively, we get:
\begin{equation}
P(t)=e^{A^Tt}{P}(0)e^{At}
\end{equation}
Note: I used the fact that $e^{-A^T\tau}$ and $A^T$ commute, i.e. $A^Te^{-A^T\tau}=e^{-A^T\tau}A^T$
A: For small time increment you get approximately
$$
P(t+dt)=P(t)+(A^TP(t)+P(t)A)dt+O(dt^2)=(I+A^T\,dt)P(t)(I+A\,dt)+O(dt^2)
$$
which means that from time step to time step you get an accumulation of these factors on both sides,
$$
P(N\,dt)=(I+A^T\,dt)^NP(t)(I+A\,dt)^N+O(Ndt^2)
$$
Now if $N\,dt=\Delta t$ we get, using $(I+B/N)^N=\exp(B)+O(B^2/N)$
$$
P(t+Δt)=e^{A^T\,Δt}P(t)e^{A\,Δt}+O(Δt\,dt)
$$
so that for $dt\to 0$ one gets exactly the claimed solution form.
A: By vectorizing $P$ and using the Kronecker product, similar to a method which can be used to solve Sylvester equations, then the differential equation can also be written as
$$
\frac{d}{dt}\,\text{vec}\,P(t) = \underbrace{\left(I \otimes A^\top + A^\top\otimes I\right)}_{M}\,\text{vec}\,P(t), \tag{1}
$$
which has the solution
$$
\text{vec}\,P(t) = e^{M\,t}\,\text{vec}\,P(0). \tag{2}
$$
By using the mixed-product property of the Kronecker product it can be shown that $I \otimes A^\top$ commutes with $A^\top \otimes I$. This commuting property allows you to write the matrix exponential as
$$
e^{M\,t} = e^{(I\,\otimes\,A^\top)\,t}\,e^{(A^\top\,\otimes\,I)\,t}. \tag{3}
$$
By using the definition of a matrix exponential and again the mixed-product property of the Kronecker product then is can also be shown that
$$
e^{X\,\otimes\,I} = e^{X} \otimes I, \quad
e^{I\,\otimes\,X} = I \otimes e^{X},
$$
thus 
$$
e^{M\,t} = \left(I \otimes e^{A^\top t}\right) \left(e^{A^\top t} \otimes I\right). \tag{4}
$$
Substituting $(4)$ into $(2)$ gives
\begin{align}
\text{vec}\,P(t) &= \left(I \otimes e^{A^\top t}\right) \left(e^{A^\top t} \otimes I\right) \text{vec}\,P(0) \\
&= \left(I \otimes e^{A^\top t}\right) \text{vec}\left(P(0)\,e^{A\,t}\right) \\
&= \text{vec}\left(e^{A^\top t}\,P(0)\,e^{A\,t}\right)
\end{align}
thus
$$
P(t) = e^{A^\top t}\,P(0)\,e^{A\,t}.
$$
