<s>[INST] How do I find all solutions of $X^2 = I_2$ in $M_2(\mathbb{N})$? I know that the solutions of the equation:
$$X^2 = I_2$$
in $M_2(\mathbb{N})$ are the $2$x$2$, natural, involutory matrices:
$$ X_1 = 
    \begin{pmatrix}
    1 & 0\\
    0 & 1\\
    \end{pmatrix}
$$
and
$$ X_2 = 
    \begin{pmatrix}
    0 & 1\\
    1 & 0\\
    \end{pmatrix}
$$

How can I arrive at these solutions by calculations?

 [/INST] Take $$X=\begin{bmatrix}a& b \\ c &d
\end{bmatrix} $$
with $a,b,c,d \in \mathbb{N}$
Now
$$X^2=\begin{bmatrix}a^2+bc & ab+bd \\ ac+cd & bc+d^2
\end{bmatrix} $$
and we have a system of equations in $\mathbb{N}$:
$$\begin{cases}
a^2+bc=1\\
ab+bd=0 \\
ac+cd=0 \\
bc+d^2=1
\end{cases}$$
For $b(a+d)0=0$, we have either $b=0$ or $a=d=0$ since we are in $\mathbb{N}$.
Suppose $b=0$, then
$$X^2=\begin{bmatrix}a^2 & 0 \\ ac+cd & d^2
\end{bmatrix} $$
It follows that $a=d=1$ and $c=0$. Thus $$X=\begin{bmatrix}1& 0 \\ 0 &1
\end{bmatrix} $$
Suppose now $a=d=0$, then
$$X^2=\begin{bmatrix}bc & 0 \\ 0 & bc
\end{bmatrix} $$
and it follows that $b=c=1$. Thus $$X=\begin{bmatrix}0& 1 \\ 1 &0
\end{bmatrix} $$
</s>