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?

 A: 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} $$
A: Hint If we write a generic element of $M_2(\Bbb N)$ as $$X := \pmatrix{a&b\\c&d} ,$$ taking the $(1, 2)$ element of both sides of $$X^2 = I_2$$ gives that $(a + d) b = 0$, so either:


*

*$a + d = 0$ (and hence $a = d = 0$), or

*$b = 0$.



Additional hint In the former case, for example, the condition simplifies to $$b c I_2 = I_2 , $$ which implies $b = c = 1$, that is, $$A = \pmatrix{0&1\\1&0} .$$

A: Clearly, this (more than elementary) question has not any interest. One wonders why our ayatollahs, so quick to close some excellent questions, have not rushed to close this one. That is more interesting -and not difficult to prove- is the following generalization
EDIT.
$\textbf{Proposition}$. Let $k$ be a positive integer. Find the matrices $X\in M_n(\mathbb{N})$ satisfying $X^k=I_n$.. 
$\textbf{Proof}$. Note that $X^{-1}=X^{k-1}\in M_n(\mathbb{N})$. More generally, when  $X,X^{-1}\in M_n(\mathbb{N})$, it is not difficult to deduce that $X$ is a permutation. Then, here, the set of solutions is the set of permutations of order $p$ where $p$ varies through the divisors of $k$. Recall that a permutation of order $p$ is a product of disjoint cycles $c_1\circ\cdots \circ c_q$ where $p=lcm(length(c_1),\cdots,length(c_q))$.
