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.