Solve a Diophantine equation with three variables How can I find all integral solutions of
$$4xy-x-y+20=z^2$$
I know two solutions are
$(x,y,z)=(-1,-3,\pm6)$. Is there a way to find new, or all integral solutions from this one?
This is not exhaustive. It is only a way given one solution, generate a family of solutions.

Let $\Lambda(x,y,z)$ be the $2\times 2$ symmetric matrices $\begin{bmatrix}4x - 1 &amp; 2z \\ 2z &amp; 4y - 1\end{bmatrix}$.  

The equation at hand can be rewritten as
$$\det \Lambda(x,y,z) = (4x-1)(4y-1)-(2z)^2 = -79\tag{*1}$$

For any $2\times 2$ matrix $P$ with integer coefficients and $\det P = \pm 1$, $P^T \Lambda(x,y,z) P$ will be a $2 \times 2$ symmetric matrix with integer coefficients. It has same determinant as $\Lambda(x,y,z)$. 
If $(x,y,z)$ is a solution of  $(*1)$ and we can find integers $x',y',z'$ such that $P^T \Lambda(x,y,z) P = \Lambda(x',y',z')$, then $(x',y',z')$ will be another solution for $(*1)$.

This prompts us to look for suitable $P$ from $SL(2,\mathbb{Z})$.
$SL(2,\mathbb{Z})$ are generated by following matrices:
$$
L = \begin{bmatrix}1 &amp; 0 \\ 1 &amp; 1\end{bmatrix},\quad
U = L^T = \begin{bmatrix}1 &amp; 1 \\ 0 &amp; 1\end{bmatrix}
\quad\text{ and }\quad
J = \begin{bmatrix}0 &amp; -1 \\ 1 &amp; 0\end{bmatrix}$$
It is easy to check $J^T \Lambda(x,y,z) J = \Lambda( y, x, -z )$, so $J$ doesn't give us anything interesting.
However, $L$ and $U$ don't disappoint us. We find

$$\begin{align}
U^{2k} \Lambda(x,y,z) L^{2k} &amp;= \Lambda(x + 2kz + k^2(4y-1), y, z + k(4y-1))\\
L^{2k} \Lambda(x,y,z) U^{2k} &amp;= \Lambda(x, y+2kz + k^2(4x-1), z + k(4x-1))
\end{align}
$$

This means if $(x,y,z)$ is an integral solution for $(*1)$, then for any integer $k$, Both $$\begin{align}
&amp; (x + 2kz + k^2(4y-1), y, z + k(4y-1))\\
\text{ and }\quad&amp;(x, y+2kz + k^2(4x-1), z + k(4x-1))
\end{align}
$$ are solutions of $(*1)$.

For example, if one start from the solution $(-1,-3,6)$, we immediately obtain
following two parametric families of solution:

$$(-1 + 12k -13k^2, -3, 6-13k)\quad\text{ and }\quad (-1, -3+12k-5k^2, 6-5k)$$

More solutions can be constructed in this manner. However, I doubt this type of
construction cover all possible solution. I hope this is at least a start.
Solving for $y$, we get
$$ y = \frac{z^2+x-20}{4x-1}$$
That is, we want $z^2 + x - 20$ to be divisible by $4x-1$.  If $m = 4x-1$, 
$z^2 + x - 20 = ((2 z)^2 + m - 79)/4$, so what you need is that $79$ is a square mod $m$, where $m \equiv -1 \mod 4$.  

The first few $m$ are 
$ 3, 7, 15, 27, 35, 39, 43, 47, 59, 63$, corresponding to $x = 1, 2, 4, 7, 9, 10, 11, 12, 15, 16$.  Then for example, for $x = 16$, $m = 63$, $2z \equiv \pm 4$ or $\pm 31 \mod 63$, making $z \equiv \pm 2$ or $\pm 16 \mod 63$.  If you take $z = 16$, $y = (16^2 + 16 - 20)/63 = 4$.

In particular, if $m$ is a prime $\equiv -1 \mod 4$, by quadratic reciprocity $79$ is a square mod $m$ if and only if $m$ is $79$ or is not a square mod $79$.  By Dirichlet's theorem we get infinitely many primes of this form.