Finding all matrices $B$ such that $AB=BA$ for a fixed matrix $A$

Here's another approach. Notice that if $A$ and $B$ commute and if $C$ is any invertible matrix, then $C^{-1}AC$ and $C^{-1}BC$ commute. This is because $$(C^{-1}AC)( C^{-1}BC) = C^{-1}ABC = C^{-1}BAC = (C^{-1}BC)(C^{-1}AC).$$

If we let $C$ be the matrix consisting of eigenvectors of $A$, then one can calculate that $$CAC^{-1} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0\\ 0 & 0& 2\end{bmatrix}.$$

So, what commutes with $CAC^{-1}$? Multiplying a matrix on the left by $CAC^{-1}$ multiplies the rows by $1$, $1$, and $2$ respectively. Multiplying a matrix on the right by $CAC^{-1}$ multiplies the columns by $1$, $1$, and $2$ respectively. The only way these can be equal is if the spots which are multiplied by both $1$ and $2$ are $0$. It follows that all of the matrices which commute with $CAC^{-1}$ have the form $$\begin{bmatrix} a & b & 0\\ c&d&0\\ 0&0&e\end{bmatrix}$$ and one easily checks that all of these do in fact commute with $C^{-1}AC$.

To turn this back into an answer for $A$ (instead of $C^{-1}AC$), multiply this general form on the left by $C$ and on the right by $C^{-1}$.

Doing this on maple, I find that the matrices which commute with $A$ are precisely those of the form $$\begin{bmatrix} a+c & \frac{1}{3}(a+c-b-d) & 0 \\ -3c & -c+d & 0 \\ -3a+3e & -a+b+e & e\end{bmatrix}. $$

Let $$B = \begin{bmatrix}b_{11} & b_{12} & b_{13}\\ b_{21} & b_{22} & b_{23}\\ b_{31} & b_{32} & b_{33} \end{bmatrix}$$ We have $$A = \begin{bmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 3 & 1 & 2 \end{bmatrix}$$ Hence, \begin{align} AB & = \begin{bmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 3 & 1 & 2 \end{bmatrix}\begin{bmatrix}b_{11} & b_{12} & b_{13}\\ b_{21} & b_{22} & b_{23}\\ b_{31} & b_{32} & b_{33} \end{bmatrix}\\ & = \begin{bmatrix}b_{11} & b_{12} & b_{13}\\ b_{21} & b_{22} & b_{23}\\ 3b_{11} + b_{21} +2b_{31} & 3b_{12} + b_{22} +2b_{32} & 3b_{13} + b_{23} +2b_{33} \end{bmatrix} \end{align} \begin{align} BA & = \begin{bmatrix}b_{11} & b_{12} & b_{13}\\ b_{21} & b_{22} & b_{23}\\ b_{31} & b_{32} & b_{33} \end{bmatrix}\begin{bmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 3 & 1 & 2 \end{bmatrix}\\ & = \begin{bmatrix}b_{11} + 3b_{13} & b_{12} + b_{13} & 2b_{13}\\ b_{21} + 3b_{23} & b_{22} + b_{23} & 2b_{23}\\ b_{31} + 3b_{33} & b_{32} + b_{33} & 2b_{33}\end{bmatrix} \end{align} Hence, we get that $$\begin{bmatrix}b_{11} & b_{12} & b_{13}\\ b_{21} & b_{22} & b_{23}\\ 3b_{11} + b_{21} +2b_{31} & 3b_{12} + b_{22} +2b_{32} & 3b_{13} + b_{23} +2b_{33} \end{bmatrix} = \begin{bmatrix}b_{11} + 3b_{13} & b_{12} + b_{13} & 2b_{13}\\ b_{21} + 3b_{23} & b_{22} + b_{23} & 2b_{23}\\ b_{31} + 3b_{33} & b_{32} + b_{33} & 2b_{33}\end{bmatrix}$$ This gives us $$b_{13} = b_{23} = 0$$ Hence, we have $$B = \begin{bmatrix} b_{11} & b_{12} & 0\\ b_{21} & b_{22} & 0 \\ 3(b_{33} - b_{11})-b_{21} & b_{33} - b_{22} - 3b_{12} & b_{33}\end{bmatrix}$$ where $b_{ij} \in \mathbb{R}$.