How to solve $AX=XB$ for $X$ Matrix? I have two symmetric $3\times 3$ matrices $A, B$. I am interested in solving the system
$$AX= XB$$
Is there a way this is usually done? The matrices are not necessarily non singular.
 A: You can use vectorization and Kronecker products.
Assume we have the equation
$$\bf AXB=C$$
$$({\bf B}^T \otimes {\bf A})\text{vec}({\bf X}) = \text{vec}({\bf AXB}) = \text{vec}({\bf C})$$
Now how can we modify it to get what we want?
A: As stated in the comment from Christoph, here is a direct approach solving $AX=XB$ given both $A$ and $B$ are square matrices of order $k$.

Proposition. Let $A=(a_{ij})_{k\times k}$, $B=(b_{ij})_{k\times k}$ and the unknown $X=(x_{ij})_{k\times k}$. We also let $$Y=(x_{11},x_{21},\cdots,x_{k1},x_{12},x_{22},\cdots,x_{k2},\cdots,x_{1k},x_{2k},\cdots,x_{kk})^{T}$$
  be the column vector with dimension $k^2$. Then, the solution of $AX=XB$ is given by the solution of the homogeneous equation $QY=0$, where $Q$ is the block matrix
  $$Q=\begin{bmatrix}
Q_{11}&Q_{12}&\cdots&Q_{1k}\\
Q_{21}&Q_{22}&\cdots&Q_{2k}\\
\cdots&\cdots&\cdots&\cdots\\
Q_{k1}&Q_{k2}&\cdots&Q_{kk}\\
\end{bmatrix}$$
  with $Q_{ii}=A-b_{ii}I$ ($i=1,2,\dots,k$ and $I$ is the identity matrix of order $k$) and $Q_{ij}=-b_{ji}I$ ($i=1,2,\dots,k$).

The proof follows by direct multiplication of matrices. We get the system of linear equations about every elements in $X$ and write the equation as the form $QY=0$. Since there always has solutions for the  homogeneous equation $QY=0$, we get all the solutions for $AX=XB$.

Example. Solve $AX=XB$, where $A=\begin{bmatrix}-4&1\\-9&2\end{bmatrix}$ and $B=\begin{bmatrix}-2&1\\-1&0\end{bmatrix}$.
From the proposition, we have
$$Q=\begin{bmatrix}
-2&1&1&0\\
-9&4&0&1\\
-1&0&-4&1\\
0&-1&-9&2\\
\end{bmatrix}\xrightarrow{\text{row operation}}
\begin{bmatrix}
1&0&4&-1\\
0&1&9&-2\\
0&0&0&0\\
0&0&0&0\\
\end{bmatrix}=P.$$
Since $PY=0$ has the same solution to $QY=0$, we get
$$Y=(-4u+v,-9u+2v,u,v)^T,\qquad \forall u,v\in\Bbb C.$$ 
Therefore,
$$X=\begin{bmatrix}
-4u+v&u\\
-9u+2v&v\\
\end{bmatrix}=
u\begin{bmatrix}-4&1\\-9&0\end{bmatrix}+v\begin{bmatrix}1&0\\2&1\end{bmatrix},\qquad \forall u,v\in\Bbb C.$$

From the example we can see the steps are


*

*express $Q$ in terms of $A$ and $B$, apply row operations to obtain $P$,

*solve $Y$ for $PY=0$,

*write down the elements in $X$ (corresponding to the elements in $Y$). 
A: This is a partial answer, but for a generalized setting where $A,B\in M_n(\Bbb{K})$ for some field $\Bbb{K}$.  Suppose that $A$ and $B$ are diagonalizable over $\Bbb K$.  Let $a_1,a_2,\ldots,a_n$ be a complete set of eigenvectors of $A$ with eigenvalues $\lambda_1,\lambda_2,\ldots,\lambda_n$, and $b_1,b_2,\ldots,b_n$ a complete set of eigenvectors of $B^t$ with eigenvealues $\mu_1,\mu_2,\ldots,\mu_n$.  Then, the matrices $C_{i,j}=a_ib_j^t$ form a basis of $M_n(\Bbb{K})$.
Consider the linear map $T:M_n(\Bbb{K})\to M_n(\Bbb{K})$ given by $T(X)=AX-XB$.  If $X=\sum_{i,j}x_{i,j}C_{i,j}$, then
$$T(X)=\sum_{i,j}x_{i,j}(\lambda_i-\mu_j)C_{i,j}.$$
In particular, $T(X)=0$ if and only if $x_{i,j}=0$ whenever $\lambda_i\neq \mu_j$.  Therefore, 
$$\ker T=\operatorname{span}\big\{C_{i,j}:\lambda_i=\mu_j\big\}.$$
In your case, the matrices $A$ and $B$ are real symmetric (you didn't specify the field, so I guess the field is the reals), and so they are diagonalizable.  Therefore, you can use my approach.
However, if the field is $\Bbb{C}$, then a symmetric matrix may not be diagonalizable.  I do not know the answer if $A$ or $B$ is not diagonalizable.  But the wiki link given by Tianlalu may help.

After I made some calculations, I found this result.  Assume that the characteristic polynomials of $A$ and $B$ split into linear factors over $\Bbb{K}$.  Suppose that $A$ has $k$ Jordan blocks of sizes $r_1,r_2,\ldots,r_k$ with eigenvalues $\lambda_1,\lambda_2,\ldots,\lambda_k$, and suppose that $B$ has $l$ Jordan blocks of sizes $s_1,s_2,\ldots,s_l$ with eigenvalues $\mu_1,\mu_2,\ldots,\mu_l$.   That is, we have the following proposition.

Proposition:  Let $T:M_n(\Bbb{K})\to M_n(\Bbb{K})$ be defined by $T(X)=AX-XB$, where $A$ and $B$ are fixed elements of $M_n(\Bbb{K})$.  Then, the linear map $T$ is diagonalizable over $\Bbb K$ if and only if both $A$ and $B$ are diagonalizable over $\Bbb{K}$.

Write $(a_i^1,a_i^2,\ldots,a_i^{r_i})$ for a generalized eigenvector sequence of $A$ in the $i$th block, that is,
$$Aa_i^p=\lambda_ia_i^p+a_i^{p+1}$$
for $p=1,2,\ldots,r_i$, with $a_i^{r_i+1}=0$.  Similarly, write $(b_j^1,b_j^2,\ldots,b_j^{s_j})$ for a generalized eigenvector sequence of $B$ in the $j$th block, that is,
$$Bb_j^q=\mu_jb_j^q+b_j^{q+1}$$
for $q=1,2,\ldots,s_j$, with $b_j^{s_j+1}=0$.  Then, for fixed $i=1,2,\ldots,k$ and $j=1,2,\ldots,l$, the span of $C_{i,j}^{p,q}=a_i^p(b_j^q)^t$ with $p$ and $q$ ranging from $\{1,2,\ldots,r_i\}$ and $\{1,2,\ldots,s_j\}$ is a direct sum of generalized eigenspaces of $T$ with eigenvalue $\lambda_i-\mu_j$.
This is because
$$\big(T-(\lambda_i-\mu_j)\big)^hC_{i,j}^{p,q}=\sum_{r=0}^h(-1)^r\binom{h}{r}C_{i,j}^{p+h-r,q+r},$$
where $C_{i,j}^{p,q}=0$ if $p>r_i$ or $q>s_j$.  Therefore, if $h>r_i-p+s_j-q$, then $$\big(T-(\lambda_i-\mu_j)\big)^hC_{i,j}^{p,q}=0.$$  In characteristic $0$, there are $m_{i,j}=\min\{r_i,s_j\}$ generalized eigenspaces and the dimensions of the generalized eigenspaces are
$$r_i+s_j-1,r_i+s_j-3,r_i+s_j-5,\ldots,r_i+s_j+1-2m_{i,j}.$$ 
This shows that
$$\operatorname{span}\{C^{r_i,s_j}_{i,j}:\lambda_i=\mu_j\big\}\subseteq\ker T\subseteq \operatorname{span}\{C_{i,j}^{p,q}:\lambda_i=\mu_j\big\}.$$
For $X\in M_n(\Bbb{K})$, we write $X=\sum_{i,j,p,q}x_{i,j}^{p,q}C_{i,j}$.  If $X\in \ker T$, then $x_{i,j}^{p,q}=0$ whenever $\lambda_i\neq \mu_j$, and so
$$T(X)=\sum_{i,j,p,q}x_{i,j}^{p,q}(C_{i,j}^{p+1,q}-C_{i,j}^{p,q+1})=\sum_{i,j}^{p,q}(x_{i,j}^{p-1,q}-x_{i,j}^{p,q-1})C_{i,j}^{p,q},$$
with $x_{i,j}^{0,q}=x_{i,j}^{p,0}=0$.  This is as far as I can go, but it tells me that, unless both $A$ and $B$ are diagonalizable, we have
$$\operatorname{span}\{C^{r_i,s_j}_{i,j}:\lambda_i=\mu_j\big\}\subsetneq \ker T\subsetneq \operatorname{span}\{C_{i,j}^{p,q}:\lambda_i=\mu_j\big\}.$$
