How to solve for the matrix $X$ in the following equation $AXB + X = CD$ How to solve for the matrix $X$ in the following equation $AXB + X = CD$?
$A$ and $B$ are full rank symmetric matrices, and there is no structure to $CD$. $CD$ just could be $C$. 
 A: Without loss of generality $CD = C$. By properties of the Kronecker product, the problem is equivalent to
$$
(B^T \otimes A) vec(X) + vec(X) = vec(C)
$$
with solution
$$vec(X) = \left( B^T \otimes A + I \right)^{-1} vec(C)
$$
assuming the inverse exists. Here $vec(A)$ is the vector obtained from the matrix $A$ by stacking its columns.
Let $A$ and $B$ have eigenvalues $\mu_i, \lambda_j$. Then the inverse exists iff $\mu_i \lambda_j \ne -1$ for all $i,j$. In particular $A$ or $B$ do not need to have full rank.
A: As whuber wrote in a comment "It's a system of linear equations--solve it as you would any system."
Here's how you could do it in Octave (or MATLAB) with YALMIP (free) plus the free solver GLPK.
>> n=5;A=randn(n);B=rand(n);C=randn(n);D=rand(n); % generate random data
>> X=sdpvar(n,n,'full') % declare X to be an n by n matrix variable
>> optimize(A*X*B+X==C*D,[],sdpsettings('solver','glpk')) % find a solution to the constraint
>> value(X) % here is the solution
   -0.6236   -2.2800   -1.7939   -0.5188   -1.0156
    0.1853    1.4420    0.9496    0.4852   -0.3270
   -0.7491   -2.6045   -1.8847   -1.2238   -0.8824
   -0.8494    1.9423   -0.1014   -1.7763    0.1937
    0.4692    0.9236    0.4308    0.9344   -0.0396

