Prove that matrix can be square of matrix with real entries

Sure. Your matrix is the matrix of a half-turn around the $x$-axis. Just take a quarter of a turn around the same axis:$$\begin{pmatrix}1&0&0\\0&0&-1\\0&1&0\end{pmatrix}\text{ or }\begin{pmatrix}1&0&0\\0&0&1\\0&-1&0\end{pmatrix}.$$

Hint. Think about the geometry of the linear transformation $T$ of space that matrix represents. Then look for a transformation $S$ such that $S^2 = T$.

I will try to present a solution "without any trial and error".

Here we have a matrix $F$ in the form $\begin{bmatrix} 1 & 0_{1 \times 2} \\ 0_{2 \times 1} & A_{2 \times 2} \end{bmatrix}$ so if $G^2=F$ then $G$ can be in the form $\begin{bmatrix} \pm 1 & 0_{1 \times 2} \\ 0_{2 \times 1} & B_{2 \times 2} \end{bmatrix}$ where $B^2=A$.

Now let's concentrate on $B^2=\begin{bmatrix} -1 & 0 \\ 0 & -1 \end{bmatrix} $.

We have two matrix equations $B^2=-I$ i.e.
$B^2+I=0$

and general equation from Cayley-Hamilton theorem for $ 2 \times 2$ matrices $B^2-\text{tr}(B)B+\det(B)I=0$.

Comparing both equations we obtain

$\text{tr(B)}=0$ , $ \det(B)=1$.

So if we denote $ B=\begin{bmatrix} a & b \\c & d\end{bmatrix} $ then $d=-a$ and consequently $-a^2-bc=1$.
These conditions are sufficient to obtain a solution not only with real numbers, but even with integer values.

For example
if $a=0$ then $b=1$, $c=-1$

($b,c$ should have always opposite signs because -$a^2-bc$ has to be positive - possible other solution $-1,1$, also $a$ and $d$ can be switched)

if $a=1$ then $b=2$, $c=-1$
if $a=2$ then $b=5$, $c=-1$
if $a=3$ then $b=5$, $c=-2$
...
if $a=8$ then $b=5$, $d=-13$ , etc... infinite number of solutions - for every integer $a$ we can find proper integer values of $b$ and $c$ from $-bc=a^2+1$..

Let's check the last listed solution.
Indeed

$ \begin{bmatrix} 8 & 5 \\-13 & -8\end{bmatrix} \begin{bmatrix} 8 & 5 \\-13 & -8\end{bmatrix} = \begin{bmatrix} 8 \cdot {8} - 5\cdot 13 & 8\cdot 5 - 5 \cdot 8 \\ -13\cdot 8 + 8 \cdot 13 & -13\cdot 5 +8\cdot 8\end{bmatrix} = \begin{bmatrix} -1 & 0 \\0 & -1\end{bmatrix} $ .

Here is a practical computational way of solving this, which has the advantage of generalizing to more difficult problems for which the "tricks' used in the other answers may not apply. This does not constitute a mathematical "proof", but it is a constructive and practical applied math solution. Solution using YALMIP under MATLAB.

x=sdpvar(3,3,'full')
optimize(x^2==[1 0 0;0 -1 0;0 0 -1],[],sdpsettings('solver','baron'))

disp(value(x))
    1.0000         0         0
         0   -0.3562   -0.9815
         0    1.1481    0.3562

disp(value(x^2))
    1.0000         0         0
         0   -1.0000         0
         0         0   -1.0000

Now to show the power of this approach, let's say we want to find a square root, such that the (3,3) element is minimized, subject to the constraint that all elements of the square root have magnitude less than or equal to 5.

x=sdpvar(3,3,'full')
optimize([x^2==[1 0 0;0 -1 0;0 0 -1],-5<=x(:)<=5],x(3,3),sdpsettings('solver','baron'))

disp(value(x))
   -1.0000         0         0
         0    4.8990   -5.0000
         0    5.0000   -4.8990

disp(value(x^2))
    1.0000         0         0
         0   -1.0000         0
         0         0   -1.0000

Think of the construction of $\mathbf C$ via $2\times2$-matrices: if $J=\begin{bmatrix}0&-1\\1&0\end{bmatrix}$, then $$J^2=-I=\begin{bmatrix}-1&0\\0&-1\end{bmatrix}.$$ Thus, using block multiplication, one obvious solution is $$B=\begin{bmatrix}1&0&0\\0&0&-1\\0&1&0\end{bmatrix}.$$