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?
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} $$