Aquí es un enfoque general. Ya que el producto de la suma de dos cuadrados es en sí la suma de dos cuadrados, entonces,
$$\tag{1}(4x^2+1)(4y^2+1) = 4z^2+1$$
es equivalente a,
$$\tag{2}(2x+2y)^2+(4xy-1)^2 = 4z^2+1$$
La solución completa a la forma,
$$\tag{3}x_1^2+x_2^2 = y_1^2+y_2^2$$
está dada por la identidad,
$$\tag{4}(ac+bd)^2 + (bc-ad)^2 = (ac-bd)^2+(bc+ad)^2$$
Entonces, uno puede equiparar los términos de (2) y (4), resolver por {x, y, z}, con {a, b, c, d} elegido tal que uno de los términos en el lado derecho es igual a la unidad.
EDITADO MUCHO MÁS TARDE:
En respuesta a sus preguntas, vamos a tener una solución más simple a (3),
$$\tag{5}(6n+2)^2+(6n^2+4n-1)^2=(6n^2+4n+2)^2+1$$
Equiparar las condiciones de (2) y (5) y nos encontramos con que,
$$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$$
Para deshacerse de la $\sqrt{N}$ y resolver de la forma,
$$an^2+bn+c^2 = \square$$
simplemente se elige,
$$n = \frac{-2cuv+bv^2}{u^2-av^2}$$
para arbitrario {u, v}. Por supuesto, como usted desea entero n, usted tiene que resolver el denominador como la ecuación de Pell $u^2-av^2 = \pm 1$.
En resumen, y después de la simplificación, un número infinito de entero de soluciones,
$$(4x^2+1)(4y^2+1) = 4z^2+1$$
está dada por la más simple,
$$x = (u-3v)(u-v)$$
$$y = 2uv$$
$$z = (u^2-2uv+3v^2)^2$$
donde,
$$u^2-3v^2=1$$
P. S. es bastante fácil de encontrar otras soluciones similares a (5), y adecuadas tendría otras ecuaciones de Pell.