I hope the answer below is somewhat helpful.

Let me first summarize some basic facts.

It is known that the equation

\begin{equation*}
 AX + XA^T = B,
\end{equation*}

has a unique solution if the matrix $A$ is *positively stable* (i.e., has spectrum in the right half plane). If $A$ is diagonal with entries $a_1,\ldots,a_n$, then the solution to the equation can be given in closed form

\begin{equation*}
 X = D \circ B,
\end{equation*}
where $D$ is a diagonal matrix with entries $1/(\bar{a}_i+a_j)$.

In the more general case, for positively stable $A$, the solution to the above equation can be represented as

\begin{equation*}
X = \int_0^\infty e^{-tA}B(e^{-tA})^Tdt
\end{equation*}

But that does not seem to be computationally that nice.

If $n$ is largish, one can still solve the linear system written using tensor products by using an iterative algorithm for solving the linear system, as long as the iterative algorithm (e.g., conjugate gradient, or other related methods) depends on just "matrix-vector" products. Because you would need to only compute $(A \otimes I + I \otimes A)x$ several times, and that can be done using matrix multiply without actually forming the tensor products.