$(x^2+1)(y^2+1)=z^2+1$ 
Prove that $$(x^2+1)(y^2+1)=z^2+1$$ has infinitely many solutions when $x,y,z \in \Bbb N$

I couldn't even find one solution for this equation.
Any ideas on how to prove it? 
Thanks in advance.
 A: $$\boxed{\color{red}{(x^2+1)(4x^4+1)= (2x^3+x)^2+1}}$$

Some motivation for that "proof without words"  $$z^2 = x^2+y^2+x^2y^2$$
Let $y=xt$ then we have $$z^2 = x^2(1+t^2+x^2t^2)$$
so if we put $t=2x$ we get $$z^2 = x^2(1+4x^2+4x^4)=x^2(2x^2+1)^2$$
So we have triples $(x,2x^2,2x^3+x)$ where $x$ is arbitrary integer.
A: Actually we can find all the solutions $(x,y,z) \in {\mathbb Z}^3$ of the given equation $(x^2+1)(y^2+1)=z^2+1$. Note right from the start that if $(x,y,z)$ is a solution, so is any $(x', y', z')$ obtained by changing independently the signs of $x,y,z$. It will be convenient to introduce the Gaussian ring $\mathbb Z[i]= \{a+bi\mid a,b \in \mathbb Z\}$, and the norm map $N:\mathbb Q(i)\to \mathbb Q$, defined by $x+yi \to (x+yi)(x-yi)=x^2+y^2$, with $x,y \in \mathbb Q$. This norm is obviously a multiplicative function.
A solution $(x,y,z)\in {\mathbb Z}^3$ of our equation is s.t. $N(x+i)N(y+i)=N(z+i)$, or equivalently, $N((x+i)(y+i)/(z+i))=1$. In a Galois extension with cyclic group generated by $\sigma$, an element $\alpha$ has norm $1$ iff it is of the form $\beta/\sigma(\beta)$ : this is Hilbert's thm.90, which can be shown "by hand" in the quadratic case. So $(x+i)(y+i)/(z+i)=(m+ni)/(m-ni)$ here, with $(0,0)\neq (m,n) \in \mathbb Q$ a priori (we could take $m,n\in \mathbb Z$ because of homogeneity). Clearing denominators and identifying the real and imaginary parts, we get a "parametrization" system consisting of two  equations, $ma+nb=0, mb-na=0$, with $a=z+xy-1$ and $b=x+y-1$. Since $(m,n)\neq (0,0)$, the determinant $a^2+b^2$ is null, i.e. $a=b=0$, so $m,n$ are arbitrary and play no role.
Summarizing : all the solutions in $\mathbb Z^3$ of our diophantine equation consist of all the triples $(x,y,z)$ with  $y=-x+1, z=-x^2+x-1$, and all the triples obtained by changing independently the signs of $x,y,z$. In particular, we recover the triples $(a,a+1,a^2+a+1)$ given by @Don Thousand. 
A: $$(n^2+1)((n+1)^2+1)=(n^2+n+1)^2+1$$
