This is probably a really basic question, but we are stuck and the usual keyword lead to the normal "Solving linear equations with matrices"/Gaussian-elimination pages ...
I have an equation $A X B + C X D = E$, with $A,B,C,D$ being square matrices and thus, dim(X)=dim(E). Now, with $B=D$ this would be trivial to solve, but in that form, I have no clear idea how to do it easily. It should be possible though, since if I multiply everything, I have as many equations as variables and everything is linear. So I could reshape the matrices to vectors, work with use sparse matrices and so on, but that doesn't sound elegant.
If it makes it somewhat easier, in my case, $B$ is the identity matrix and D only consists of ones. In the end, I want to solve this using MATLAB.
Let $\def\vec{\mathop{\rm vec}}\vec\colon \def\Mat{\mathop{\rm Mat}}\Mat(n)\to \mathbb R^{n^2}$ denote vectorisation (columnwise, that is $A_{ij} = \vec(A)_{i+n(j-1)}$) and $\otimes$ the Kronecker product of matrices (
kronin matlab), then we have $$ \vec(AXB) = (B^t \otimes A) \vec(X)$$ Applying this to your equation gives $$ \vec(E) = \vec(AXB) + \vec(CXD) = (B^t \otimes A + C^t\otimes D)\vec(X) $$ Now solve for $\vec(X)$ an reshape.