Solving a three-variable Diophantine

[SeventhSigma]

I have the following equation

$$(4x^2+1)(4y^2+1) = (4z^2+1)$$

For positive, nonzero integers x and y (and thus z). I am having difficulty figuring out a good method/algorithm for calculating solutions to this equation. Any thoughts?

 

[RamaWolf]

One nice solution: x=56, y=209, z=23409

 

[RamaWolf]

A 'good deal' of the solutions are caught by:

let x be a natural number and y = 4 x^{2},

then we have z = x (2 y +1)^{2}

for example (x,y,z) = (1,4,9), (2,16,66), (3,36,219), ..., (17,1156,39321)