<s>[INST] Find all non-similar solutions of matrix equation 
Find all unique $($not conjugate to each other by an element of $GL_2(\mathbb{Z}))$ matrices $A \in M_2(\mathbb{Z})$, such that $A ^ 2 - 4A - I = 0$, where $I$ is the identity matrix.

From the equation it can be easily seen that matrix is a solution iff it has trace equal to $4$ and determinant equal to $-1$. The main problem is to find all conjugacy classes. Since $\mathbb{Z}$ is not a field, Jordan, Frobenius and Smith normal forms are not applicable. I also tried to brute-force the problem and solve the equation $CA = BC$, where $A$ and $B$ are two distinct solutions and $C \in GL_2(\mathbb{Z})$, but this boils down to the system of 3 diophantine equations with 4 variables, one of which is not even linear (as in this paper, page 3) and I have no idea how to check existence of solutions for that.
Another approach that I tried (again, without luck) was to decompose $A$ and $B$ into a product of $GL_2(\mathbb{Z})$ generators and check that they differ by a cyclic permutation (as listed here, point 6) but I was not able to do that in general for this family of matrices.
I have also found a theorem (in this presentation by Svetlana Katok, slide 10), which is very cryptic to me, but I think it is impossible to use in practice.
I feel that the answer is $\pmatrix{1 & 2 \\ 2 & 3}$ and $\pmatrix{0 & 1 \\ 1 & 4}$ because I could not find any counterexamples.
Is there any simple way to find all conjugacy classes?
 [/INST] Your guess is correct. Here is an elementary proof. For $b\neq 0$, let
$$
F(a,b)=\begin{pmatrix} a && b \\ -\frac{a^2-4a-1}{b} && 4-a \end{pmatrix} \tag{1}
$$
Then the initial matrix $A$ is necessarily of the form $F(a,b)$ for some integers $a$ and $b$ (such that $b$ divides $a^2-4a-1$). All the conjugates of $A$ will also be of this form.
The first step is then to find, when $b$ is large enough, a suitable conjugate $F(a',b')$ of $F(a,b)$ with $|b'| \lt |b|$ ; this will allow us to argue by induction on $|b|$.
Let $t\in{\mathbb Z}$ and
$$ P=\begin{pmatrix} 1 && 1 \\ -(t+1) && -t \end{pmatrix} \in SL_2({\mathbb Z})\tag{2} $$.
Then we have (where $\ldots$ represent terms that we need not compute)
$$ 
\begin{array}{lcl}
P^{-1}AP &=& \begin{pmatrix} -t && -1 \\ \ldots && \ldots \end{pmatrix} 
        \begin{pmatrix} a && b \\ -\frac{a^2-4a-1}{b} && 4-a \end{pmatrix}
        \begin{pmatrix} 1 && 1 \\ -(t+1) && -t \end{pmatrix}  \\
        &=& \begin{pmatrix} -t && -1 \\ \ldots && \ldots \end{pmatrix} 
        \begin{pmatrix} \ldots && a-bt \\ \ldots && -\frac{a^2-4a-1}{b}+(a-4)t \end{pmatrix} \\
        &=& 
        \begin{pmatrix} \ldots && -t(a-bt)+\frac{a^2-4a-1}{b}-(a-4)t \\ \ldots && \ldots \end{pmatrix} \\
\end{array}\tag{3}$$
We have thus shown that $P^{-1}F(a,b)P$ is of the form $F(a',b')$ where
$$ 
\begin{array}{lcl}
b' &=& -t(a-bt)+\frac{a^2-4a-1}{b}-(a-4)t \\
   &=& \frac{1}{b}\bigg(-bt(a-bt)+(a^2-4a-1)-(a-4)bt\bigg) \\ 
   &=& \frac{1}{b}\bigg(-bt(a-bt)+(a^2-4a)-(a-4)bt-1\bigg) \\ 
   &=& \frac{1}{b}\bigg(-bt(a-bt)+(a-4)(a-bt)-1\bigg) \\ 
   &=& \frac{1}{b}\bigg((a-bt-4)(a-bt)-1\bigg) \\ 
   &=& \frac{1}{b}\bigg((a-bt-2)^2-5\bigg) \\ 
\end{array}\tag{4}$$
We can certainly find an integer $t$ such that $|a-2-bt|\leq \frac{|b|}{2}$, and
we then deduce from (4) that $|b'|\leq \frac{\max\bigg(5,\big|\frac{b^2}{4}-5\big|\bigg)}{|b|}$. When
$|b| \geq 7$, we have  $\max(5,\big|\frac{b^2}{4}-5\big|) \leq b^2$ and hence $|b'| \leq |b|$.
If we look a little closer, we realize that we can extend this argument : for $|b|\geq 3$ we can still find a $t$ making $|b'|\lt |b|$. Indeed, notice first that since $a^2-4a-1=0$ has no integral solutions modulo $3$, $|b|$ cannot equal $3$ or $6$. So we are left with $|b|=4$ or $5$.
When $|b|=4$, $a$ must be odd ; we can write $a=4s\pm 1$ for some integer $s$.
If $a=4s-1$, then $a-4(s-1)-2=1$ ; taking $t=sign(b)(s-1)$ therefore yields $b'=-1$ in (4).
If $a=4s+1$, then $a-4(s-1)-2=3$ ; taking $t=sign(b)(s-1)$ therefore yields $b'=1$ in (4).
When $|b|=5$, $a$ must be of the form $5s+2$ for some integer $s$, so that $a-5(s-1)-2=5$ ; taking $t=sign(b)(s-1)$ therefore yields $b'=4$ in (4).
We are now left with the case $|b|=1$ or $2$, which is the base case of the induction, and here we apply a more direct argument :
When $|b|=1$, so that $b=\varepsilon$ with $\varepsilon \in \lbrace \pm 1\rbrace$, we can write
$$
F(a,b)=Q\begin{pmatrix} 0 && 1 \\ 1 && 4 \end{pmatrix}Q^{-1}, \ \textrm{where} \ Q=\begin{pmatrix} 0 && 1 \\ \varepsilon && \varepsilon(4-a) \end{pmatrix}. \tag{5}
$$
When $|b|=2$, so that $b=2\varepsilon$ with $\varepsilon \in \lbrace \pm 1\rbrace$, we can write
$$
F(a,b)=Q\begin{pmatrix} 1 && 2 \\ 2 && 3 \end{pmatrix}Q^{-1}, \ \textrm{where} \ Q=\begin{pmatrix} 0 && 1 \\ \varepsilon && \varepsilon\frac{3-a}{2} \end{pmatrix}. \tag{6}
$$
This finishes the proof.
</s>