Solution of a linear matrix equation i'm trying to rewrite a linear matrix equation 
$$
X=A_1XA_2+A_3X^TA_4+B
$$
where $ 
X\in\mathbb{R}^{m\times n}
$ and 
$
 A_1\text{, }A_2\text{, }A_3\text{, }A_4\text{, }B  
$ are matrices of appropriate sizes.
So far i got to express the equation in terms of the components of $X$:
$$
X_{ij}-\sum_{k=1}^{m}\sum_{l=1}^{n}\left(A_{1,ik}\cdot A_{2,lj}+A_{3,il}\cdot A_{4,kj}\right)\cdot X_{kl}=B_{ij}
$$
but is there any neat way to express  this linear system in terms of a coefficient matrix $C$ and unfolded vectors $x$ and $b$
$$
Cx=b
$$
so that i can feed the problem to Matlab? I'd rather avoid for loops, since the matrices are quite big and the linear problem is framed inside a bigger, iterative computation.
Thanks!
 A: As the comment states, we can use vectorization.  In Matlab notation, $\operatorname{vec}(X)$ just is $X(:)$, and we can go from $\operatorname{vec}(X)$ to $X$ using the reshape command.  
Vectorizing both sides yields
$$
X - A_1 XA_2 - A_3 X^TA_4 = B \implies\\
\operatorname{vec}[X - A_1 XA_2 - A_3 X^TA_4] = \operatorname{vec}(B) \implies\\
\operatorname{vec}(X) - A_2^T \otimes A_1 \operatorname{vec}(X) - A_4^T \otimes A_3\operatorname{vec}(X^T) = \operatorname{vec}(B).
$$
Here, $\otimes$ denotes the Kronecker product, implemented in Matlab with kron. This is almost in the form that we want, but we need to deal with $X^T$.  For that, we need the fact that $\operatorname{vec}(X^T) = P\operatorname{vec}(X)$, where 
$$
P = \sum_{i=1}^m \sum_{j=1}^n (e_j \otimes e_i)(e_i \otimes e_j)^T = 
\sum_{i=1}^n \sum_{j=1}^m \operatorname{vec}(E_{ij})[\operatorname{vec}(E_{ji})]^T.
$$
$P$ can be built with a quick for-loop and only needs to be put together once (assuming $m,n$ are constant).  I don't know of a nicer approch here, unfortunately.
From there, we can finally write
$$
\operatorname{vec}(X) - A_2^T \otimes A_1 \operatorname{vec}(X) - (A_4^T \otimes A_3)P\operatorname{vec}(X) = \operatorname{vec}(B) \implies\\
(I - A_2^T \otimes A_1 - (A_4^T \otimes A_3)P)\operatorname{vec}(X) = \operatorname{vec}(B) \implies\\
(I - A_2^T \otimes A_1 - (A_4^T \otimes A_3)P)x = b
$$
where $x = \operatorname{vec}(X),b = \operatorname{vec}(B)$.
