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}. $$
The Cayley-Hamilton theorem states that every matrix $A$ satisfies its own characteristic polynomial; that is the polynomial for which the roots are the eigenvalues of the matrix:
$p(\lambda)=\det[A-\lambda\mathbb{I}]$.
If you view the polynomial:
$a^2-5a+6=0$,
as a characteristic polynomial with roots $a=2,3$, then any matrix with eigenvalues that are any combination of 2 or 3 will satisfy the matrix polynomial:
$A^2-5A+6\mathbb{I}=0$,
that is any matrix similar to:
$\begin{pmatrix}3 & 0\\ 0 & 3\end{pmatrix}$,$\begin{pmatrix}2 & 0\\ 0 & 2\end{pmatrix}$,$\begin{pmatrix}2 & 0\\ 0 & 3\end{pmatrix}$. Note:$\begin{pmatrix}3 & 0\\ 0 & 2\end{pmatrix}$ is similar to $\begin{pmatrix}2 & 0\\ 0 & 3\end{pmatrix}$.
To see why this is true, imagine $A$ is diagonalized by some matrix $S$ to give a diagonal matrix $D$ containing the eigenvalues $D_{i,i}=e_i$, $i=1..n$, that is:
$A=SDS^{-1}$, $SS^{-1}=\mathbb{I}$.
This implies:
$A^2-5A+6\mathbb{I}=0$,
$SDS^{-1}SDS^{-1}-5SDS^{-1}+6\mathbb{I}=0$,
$S^{-1}\left(SD^2S^{-1}-5SDS^{-1}+6\mathbb{I}\right)S=0$,
$D^2-5D+6\mathbb{I}=0$,
and because $D$ is diagonal, for this to hold each diagonal entry of $D$ must satisfy this polynomial:
$D_{i,i}^2-5D_{i,i}+6=0$,
but the diagonal entries are the eigenvalues of $A$ and thus it follows that the polynomial is satisfied by $A$ iff the polynomial is satisfied by the eigenvalues of $A$.
Best Answer
Here is what I tried. Let $A=\left(\begin{array}{cc}1&2\\3&7\\\end{array}\right)$. Then $\chi_A(x)=x^2-8x+1$ which has two roots $\lambda_1=4+\sqrt{15}$ and $\lambda_2=4-\sqrt{15}$. So there exists $P$ invertible such that $$P^{-1}AP= \left(\begin{array}{cc}\lambda_1&0\\0&\lambda_2\\\end{array}\right)=:D.$$ Now for $X$ a given matrix, $X^2=A\iff Y^2=D$ where $Y=P^{-1}XP$. Hence we just have to find $Y$ satisfying $Y^2=D$, and the only possible solutions to this are the four matrix of the form $$\left(\begin{array}{cc}\pm\sqrt{\lambda_1}&0\\0&\pm\sqrt{\lambda_2}\\\end{array}\right).$$ To see this, you can see the link given in the comments. One way to prove this would be as follows. Let $\mu_1$ and $\mu_2$ be two complex eigenvalues of $Y$ (maybe equal). Then $$Y\sim\left(\begin{array}{cc}\mu_1&\alpha\\ 0&\mu_2\\ \end{array}\right)$$ for some $\alpha$ so $$Y^2\sim\left(\begin{array}{cc}\mu_1^2&\alpha^\prime\\ 0&\mu_2^2\\ \end{array}\right).$$ From here we deduce that $\mu_1\not=\mu_2$, and for example $\mu_i^2=\lambda_i$. Futhermore, $\mu_1$ can't be a complex number (not real I mean) because we would have $\mu_2 = \bar{\mu_1}$ and $(\mu_1\mu_2)^2=det(D)=1=\vert \mu_1\vert^2$ which is not the case ($\vert \lambda_i\vert\not=1$). This shows that $Y$ has two real eigenvalues, so is diagonalizable. Because $D$ and $Y$ commute they can be simultaneously diagonalizable, which means that $Y$ is in fact diagonal, and so of the form $\left(\begin{array}{cc}\pm\sqrt{\lambda_1}&0\\0&\pm\sqrt{\lambda_2}\\\end{array}\right).$
Finally the solutions of the problem are the matrix of the form $$P\left(\begin{array}{cc}\pm\sqrt{\lambda_1}&0\\0&\pm\sqrt{\lambda_2}\\\end{array}\right)P^{-1}.$$ You just have to find $P$ (you have $P= (X_1~X_2)$ where $PX_i=\lambda_i X_i,~i=1,2$, so there are just two linear systems to solve).
I hope I haven't done any mistake, sorry for my bad English I'm not so good at it!