Since $f$ and $g$ are polynomial, using SymPy's solve_poly_system
:
>>> from sympy import *
>>> x, y, z, mu = symbols('x y z mu', real=True)
>>> L = x**2 + y**2 + z**2 + mu * (x*y - z + 1)
>>> solve_poly_system([diff(L,x), diff(L,y), diff(L,z), diff(L,mu)], x, y, z, mu)
[(0, 0, 1, 2), (-sqrt(2)*I, -sqrt(2)*I, -1, -2), (sqrt(2)*I, sqrt(2)*I, -1, -2)]
Hence, the only real solution is $(x,y,z,\mu) = (0, 0, 1, 2)$. Not very insightful, however.
Let $\mathcal L$ be the Lagrangian. Computing $\partial_x \mathcal L$, $\partial_y \mathcal L$ and $\partial_z \mathcal L$ and finding where they vanish,
$$\begin{bmatrix} 2 & \mu & 0\\ \mu & 2 & 0\\ 0 & 0 & 2\end{bmatrix} \begin{bmatrix} x\\ y\\ z\end{bmatrix} = \begin{bmatrix} 0\\ 0\\ \mu\end{bmatrix}$$
Note that the matrix is singular when $\mu = \pm 2$. Hence, we have three cases to consider.
$\color{blue}{\boxed{\mu = 2}}$
The solution set is the line parameterized by
$$\begin{bmatrix} x\\ y\\ z\end{bmatrix} = \begin{bmatrix} t\\-t\\ 1\end{bmatrix}$$
and, since, $xy - z + 1 = 0$, we obtain $t = 0$ and $\color{blue}{(x,y,z) = (0,0,1)}$.
$\color{blue}{\boxed{\mu =-2}}$
The solution set is the line parameterized by
$$\begin{bmatrix} x\\ y\\ z\end{bmatrix} = \begin{bmatrix} t\\ t\\ -1\end{bmatrix}$$
and, since, $xy - z + 1 = 0$, we obtain the equation $t^2 = -2$, which has no solution over the reals.
$\color{blue}{\boxed{\mu \neq \pm2}}$
The solution set is the line parameterized by
$$\begin{bmatrix} x\\ y\\ z\end{bmatrix} = \begin{bmatrix} 0\\ 0\\ \frac{\mu}{2}\end{bmatrix}$$
and, since, $xy - z + 1 = 0$, we obtain $\mu = 2$, which contradicts the assumption.

