There are infinitely many integer solutions, so you need to specify what you mean by "can one find all". If you mean "Are all solutions given by a finite set of easily computed parametrized families?" then the answer is yes, essentially following John R Ramsden's answer (with some minor corrections).

If we set $X = 2x+1$, $Y = 2y+1$ and $Z = 2z+1$, we get the equation $X^2 + 1 = Y^2 + Z^2$. Euler showed that all solutions have the form $(ac+bd)^2 + 1 = (ac-bd)^2 + (ad+bc)^2$ where $a,b,c,d$ range over integers satisfying $ad-bc = 1$. This set is easy to parametrize, because $a,b,c,d$ describe entries of elements $\binom{ab}{cd}$ in $SL_2(\mathbb{Z})$. We are seeking odd solutions, and this condition is equivalent to $ac+bd$ odd, which can be determined by reducing matrices mod 2.

We find that any element of $SL_2(\mathbb{Z})$ congruent to $\binom{11}{01}$, $\binom{01}{11}$, $\binom{10}{11}$, or $\binom{11}{10}$ modulo 2 yields a solution to your equation. These can be easily generated, either by using the Euclidean algorithm to generate suitable elements of $SL_2(\mathbb{Z})$ from an initial choice of coprime integers, or by multiplying an initial solution by powers of the matrices $\binom{12}{01}$ and $\binom{10}{21}$. The kernel $\Gamma(2)$ of reduction mod 2 in $SL_2(\mathbb{Z})$ becomes a free group generated by those two matrices, once we quotient by $\pm \binom{10}{01}$. This sign ambiguity nicely cancels the invariance under sign change in Euler's parametrization.