Consider the following equation with integral, nonzero $x,y,z$
$$(4x^2+1)(4y^2+1) = (4z^2+1)$$
What are some general strategies to find solutions to this Diophantine?
If it helps, this can also be rewritten as $z^2 = x^2(4y^2+1) + y^2$
I've already looked at On the equation $(a^2+1)(b^2+1)=c^2+1$
On
Here is one general approach. Since the product of the sum of two squares is itself the sum of two squares, then,
$$\tag{1}(4x^2+1)(4y^2+1) = 4z^2+1$$
is equivalent to,
$$\tag{2}(2x+2y)^2+(4xy-1)^2 = 4z^2+1$$
The complete solution to the form,
$$\tag{3}x_1^2+x_2^2 = y_1^2+y_2^2$$
is given by the identity,
$$\tag{4}(ac+bd)^2 + (bc-ad)^2 = (ac-bd)^2+(bc+ad)^2$$
One can then equate the terms of (2) and (4), solve for {x, y, z}, with {a, b, c, d} chosen such that one term on the RHS is equal to unity.
EDITED MUCH LATER:
In response to your questions, let's have a simpler solution to (3) as,
$$\tag{5}(6n+2)^2+(6n^2+4n-1)^2=(6n^2+4n+2)^2+1$$
Equate the terms of (2) and (5) and we find that,
$$x = \frac{1}{2}\big(1+3n-\sqrt{3n^2+2n+1}\big)$$
$$y = \frac{1}{2}\big(1+3n+\sqrt{3n^2+2n+1}\big)$$
$$z = (6n^2+4n+2)/2$$
To get rid of the $\sqrt{N}$ and solve the form,
$$an^2+bn+c^2 = \square$$
one simply chooses,
$$n = \frac{-2cuv+bv^2}{u^2-av^2}$$
for arbitrary {u, v}. Of course, since you want integer n, you have to solve the denominator as the Pell equation $u^2-av^2 = \pm 1$.
In summary, and after simplification, an infinite number of integer solutions to,
$$(4x^2+1)(4y^2+1) = 4z^2+1$$
is given by the rather simple,
$$x = (u-3v)(u-v)$$
$$y = 2uv$$
$$z = (u^2-2uv+3v^2)^2$$
where,
$$u^2-3v^2=1$$
P.S. It is quite easy to find other solutions similar to (5), and appropriate ones would need other Pell equations.
Let $a$ be a positive integer.
Then
\begin{align} (4a^2+1)(4((2a)^2)^2+1) &= 256a^6 + 64a^2 + 4a^2 + 1 \\ & = 4(64a^6 + 16a^4 + a^2) + 1 \\ &= 4(a^2(8a^2+1)^2)+1 \\ &= 4((8a^2+1)a)^2+1 \end{align}
so $(a, (2a)^2, (8a^2+1)a)$ is always a solution.
There are others as well.