Plug in $x = \epsilon$, $y = 3\epsilon^2$ and you will get that the function is greater than $0$ for all $\epsilon > 0$. Plug in $x = (3/4) \epsilon$, $y = \epsilon ^2$ and you will get that the function is less than $0$ for all $\epsilon > 0$. This implies that $(0,0)$ must be a saddle point because you take take $\epsilon > 0$ arbitrarily small.
After setting all three equations to $0$, multiply the first equation by $z$, the second by $x$, and the third by $xy$ to get:
\begin{align}
2xz + 2xyz &= 0 \quad (1) \\
x^3 + 2xyz &= 0 \quad (2) \\
xy^3 + 2xyz - 4xy &= 0 \quad (3) \\
\end{align}
\begin{align}
(1) - (2) &\Rightarrow 2xz - x^3 = 0 &(4) \\
(2) - (3) &\Rightarrow x^3 - xy^3 -4xy = 0 &(5) \\
(1) - (3) &\Rightarrow 2xz - xy^3 -4xy = 0 &(6) \\
\end{align}
Subtracting $(6)$ from $(5)$ gives $x^3 - 2xz = x(x^2-2z) = 0$ so that $x = 0$ or $z = x^2/2$.
In the first case, we get $2yz=0$ from the second equation, so either $y=0$ or $z=0$.
If $y=0$, then $2z-4=0$ so that $z=2$, which gives the solution $(0,0,2)$.
If $z=0$ then $y^2 = 4$ so that $y=\pm2$, which gives the solutions $(0,2,0)$, and $(0,-2,0)$.
In the second case, we get $x\neq0$ so that $z=x^2$.
Then, first original equation gives $2x(1+y) = 0$. But $x \neq 0$, so we must have $y=-1$.
Plug this into the third original equation to get $(-1)^2 + 2z - 4 = 0$, i.e. $z = \frac{3}{2}$.
Now plug $y=-1$ and $z=3/2$ into the second original equation to get $x^2 +2(-1)(\frac{3}{2}) = 0$ which gives $x^2 = 3$ or $x = \pm\sqrt 3$.
This gives the final solutions, $(-\sqrt 3, -1, \frac{3}{2})$, and $(\sqrt 3, -1, \frac{3}{2})$.
So there are three solutions: $(0, 0, 2)$, $(0,-2,0)$, $(0, 2, 0)$, $(-\sqrt 3, -1, \frac{3}{2})$, and $(\sqrt 3, -1, \frac{3}{2})$. You can verify that these work by plugging them in to the original equations.
Best Answer
Method 1 (with the Hessian): the quadratic form is
$D^2f_X(H,H)=4tr(H^TH(X^TX-I_m)+H^TXH^TX+H^TXX^TH)$
where $H\in M_{n,m}$.
When $X=0$, $D^2f_0(H,H)=-4tr(H^TH)$, that is $<0$ when $H\not=0$. Then our quadratic form is $<0$ and, since $X=0$ is a critical point of $f$, it is also a local maximum of $f$.
Method 2. Let $spectrum(X^TX)=(\sigma_i^2)$ (the singular values of $X$). Then
$f(X)=tr((X^TX-I)^2)=\sum_{i\leq m}(\sigma_i^2-1)^2$. When $X$ is in a neighborhood of $0$, then the $\sigma_i^2$ are small and
$f(X)\approx \sum_{i\leq m}(1-2\sigma_i^2)=m-2\sum_i\sigma_i^2\leq m=f(0)$ and we are done.
Remark. Of course, there are other critical points than $X=0$ and $X$ pseudo orthogonal.