Solve in integers: $y(x^2+1)=z^2+1$ Find all integer solutions to the equation
$$
y(x^2+1)=z^2+1.
$$
There is, for example, an infinite family of solutions $x=u$, $y=(uv\pm1)^2+v^2$, $z=(u^2+1)v \pm u$, $u,v \in {\mathbb Z}$, but there are also solutions outside of this family, e.g. $(x,y,z)=(8,5,18)$ or $(x,y,z)=(12,2,17)$. The question is to describe all integer solutions. Any reasonable description is ok. An algorithm generating all solutions is also ok, provided that it does not involve any search by trial and error (otherwise there is a trivial algorithm that tries all triples $(x,y,z)$ in some order). Using parametric expressions, recurrence relations, or something like ``start with this solution and apply these operations in any order'' (like generating Markov numbers via Markov tree) would be ideal, but more complicated algorithms are also possible. For example, for equation $yz=x^3+1$, there is an obvious algorithm "let $x$ be an  arbitrary integer, let $y$ be any divisor of $x^3+1$, and then let $z=(x^3+1)/y$", which I think is acceptable. The equation in question is one of the smallest/simplest ones for which I do not see any reasonable method/algorithm to describe all solutions, hence the question.
Remark: If we define, for any integer $u$, set $S(u)$ as a set of integers $0\leq r < u^2+1$ such that $\frac{r^2+1}{u^2+1}$ is an integer, then all the solutions to the equation are in the form $x=u$, $z=(u^2+1)v \pm r$ and $y=(u^2+1)v^2 \pm 2vr + \frac{r^2+1}{u^2+1}$ for $u,v \in {\mathbb Z}$ and $r \in S(u)$, but we need trial and error to construct set $S(u)$, so I do not think this is an acceptable answer.
 A: The equation says that $z^2 + 1 \equiv 0 \mod (x^2+1)$.  For each positive integer $x$, you can enumerate the square roots  of $-1$ in the integers mod $x^2+1$ as you remark in your last paragraph.  This does not require "trial and error": if you can factor $x^2+1$, you can use the Tonnelli-Shanks algorithm
to find the square roots mod each prime power in the factorization, then put them together using the Chinese Remainder Theorem.
A: This answer (which I made Community Wiki) attempts an explicit summary of how to ``list" explicitly the solutions of the given equation, given the answer of Robert Israel and the comment of David Speyer which followed it.
It's convenient to deal with the prime $2$ first, but this is easy, because the only condition imposed on $z$ to get the right power of $2$ dividing $z^{2}+1$ is that if $x$ is odd, we need $z$ odd.
If $p$ is an odd prime such that $p^{n}$ is the exact power of $p$ dividing $x^{2}+1$, then we just require that $p^{n}$ divides $z^{2} -x^{2}$, so we need $ z \equiv \pm x$ (mod $p^{n}$)  since $p$ is odd.
Notice then that, in all cases, there are $2^{k}$ possible congruences for $z$ (mod $x^{2}+1$), where $k$ is the number of distinct odd primes dividing $x^{2}+1$.
To be more explicit, if $x$ is even, let the prime factorization of $x^{2}+1$ be $x^{2}+1 = \prod_{i=1}^{k}p_{i}^{n_{i}}$, where the $p_{i}$ are distinct primes, and if $x$ is odd, let the prime factorization of $\frac{x^{2}+1}{2}$ be $\frac{x^{2}+1}{2} = \prod_{i=1}^{k}p_{i}^{n_{i}}$.
For any ordered $k$-tuple of signs $(\epsilon_{1}, \ldots, \epsilon_{k})$, solve (simultaneously) $z \equiv \epsilon_{i} x$ (mod $p_{i}^{n_{i}}$) for each $i$. If $x$ is odd, also impose the condition that $z$ is odd.
In either case, this gives $2^{k}$ distinct possibilities for the congruence of $z$ (mod $x^{2}+1$) which yield $x^{2} + 1 | z^{2}+1$, and there are no other valid congruences (where, as noted earlier, $k$ is the number of distinct odd prime divisors of $x^{2}+1$).
For each integer $r$ and each integer $m \neq 0$ let $S(m,r)$ be the set of integers $0\leq w<|m|$ such that $w^2 \equiv r \,(\text{mod } m)$. We have just explicitly described the set $S(x^2+1,-1)$ for any integer $x$. Using it, the complete solution set to our equation is given by $x=u$, $z=(u^2+1)v + r$ and $y=(u^2+1)v^2 + 2vr + \frac{r^2+1}{u^2+1}$ for $u,v \in {\mathbb Z}$ and $r \in S(u^2+1,-1)$.
A: There are two possible directions of solutions. The first one was mentioned in the previous post. It allows you to find all solutions, but it does not imply finding parametrization of solutions. This is formally an equation of two variables. And parametrization is written through solutions of the Pell equations. The difficulty of finding all solutions is that it will be necessary to consider all possible combinations and equivalent forms of the Pell equations.
A similar equation was considered there. Equation of the form.  $aX^2+bX=cY^2+dY$
Sum of two consecutive squares equals difference of two consecutive cubes
https://artofproblemsolving.com/community/c3046h1049910___4
One of the forms of solutions may look like this. $$(a^2+1)(b^2+1)=c^2+1$$
https://artofproblemsolving.com/community/c3046h1162301_almost_binary_equation
It can be reduced directly to the Pell equation. $Z^2-dR^2=c=ab$
https://artofproblemsolving.com/community/c3046h1056199_general_pell_equation
You can also consider using a general formula.   $aX^2+bXY+cY^2=f$
https://artofproblemsolving.com/community/c3046h1048219___2
$y(x^2+1)=z^2+1$
$z^2-yx^2=y-1$
Usually a replacement of this type is used. $z=Z+ax$
$Z^2+2aZx+(a^2-y)x^2=y-1$
Using a general formula, you can write a solution if the root is whole. For example so.
$$\sqrt{\frac{y-1}{1+2a+a^2-y}}=\sqrt{\frac{y-1}{(a+1)^2-y}}=1$$
$y=\frac{(a+1)^2+1}{2}$
This means that given values of $y$, there will be infinitely different values of $x,z$ that will be solutions.
Of course, there may be an option when a different representation allows you to find a simpler formula, but this approach makes it easier to solve the equation if the coefficients and their shape change. And there is no need to find a solution for an equation of another form.
