For a small project I am working on, I wish to find the solutions for $$x^2+y^2=z(4z+1)$$ in natural numbers $x,y,z$. I wish to automate finding solutions for $z$ up to a maximum value as efficient as possible, running through $z$ values from 1 and then try to find possible $x$ and $y$.
One thing that I found is that I can disregard $z \equiv 3,6,7 \mod 8$, as the sum of two squares can only be $\equiv 0,1,2,4,5 \mod 8$. But I wonder which other criteria I can use to exclude $z$ values. I also wonder whether for given $z$, which $x$ I can ignore, so that I only test those $x$ values that could give a solution.