Computing $(A\otimes I + I \otimes A)^{-1} \text{vec}B$

Your equation is equivalent to solving the following equation for $P$

$$ A\,P + P\,A^\top = B. \tag{1} $$

Such equation also solves for the infinite time (controllability) Gramian, which can also be calculated with

$$ P = -\int_0^\infty e^{A\,t}\,B\,e^{A^\top t}\,dt. \tag{2} $$

This integral should also be of size $n \times n$. As time goes to infinity $e^{A\,t}$ is not well defined when $A$ is positive semi-definite. It can also be noted that in order to be able to solve for $P$ the spectra of $A$ and $-A$ need to be disjoint. In order words $A$ can't have an eigenvalue of zero, which would imply that $A$ should be positive definite instead of positive semi-definite. When multiplying both sides of $(1)$ by minus ones keeps it equivalent to the original problem. This is also equivalent to multiplying both $A$ and $B$ by minus one, which transforms $(2)$ into

$$ P = \int_0^\infty e^{-A\,t}\,B\,e^{-A^\top t}\,dt. \tag{3} $$

Now the term $e^{-A\,t}$ does stay bounded and converges to zero as time goes to infinity. The rate of how fast this integral converges should be roughly inversely proportional to the smallest eigenvalue of $A$, however any finite time numerical approximation of this integral can be used as a lower bound for $P$.