Added: a careful proof that all solutions arise in the way described below begins with the observation that a solution $(x,y),$ with $|x|+|y|$ large, has either $x \approx (\sqrt 2 - 1) y$ or $y \approx (1-\sqrt 2)x.$ In the first case, we find $(2x-y) \approx -(3 - \sqrt 8) y.$ In the second case, $x+2y \approx (3 - \sqrt 8) x.$ Note that $(3 - \sqrt 8) \approx 0.17.$ That is, any solution can be moved until, say, $|x|+|y| \leq 10$ and all solutions in that diamond can be found.
The method for linking up all integer solutions is usually called Vieta Jumping on this site. In the infamous contest questions where this was first used, the quadratic form was of the type $x^2 - kxy + y^2.$ As it happens, there is no difficulty working the same technique when the quadratic part is $x^2 \pm kxy - y^2.$ This is a special case of the automorphism group of a (binary) quadratic form. As the coefficients of $x^2$ and $y^2$ both have absolute value $1,$ it becomes unnecessary to complete any squares or explicitly consider any Pell equations. The curious might borrow Binary Quadratic Forms by Buchmann and Vollmer. My favorite is Binary Quadratic Forms by Duncan A. Buell.
All the integer points can be produced by alternating two mappings beginning at the origin. The result is a double spiral, counterclockwise going outward, or clockwise going inward back to the origin.
Given an integer point on the spiral $(x,y),$ the neighboring spiral point with the same $x$ value, joined by a vertical line segment, is
$$ (x, 1+2x-y) $$
Given an integer point $(x,y),$ the neighboring spiral point with the same $y$ value, joined by a horizontal line segment, is
$$ (1-x-2y, y) $$
If you use on of these mapping twice in a row, you go back to where you started; this is why we alternate.
int x,y;
x = 1189; y = 2871;
cout << " $$ ( " << x << ", " << y << " ) $$ " << endl;
while ( abs(x) < 10000 && abs(y) < 10000)
{
y = 1 + 2 * x - y;
cout << " $$ ( " << x << ", " << y << " ) $$ " << endl;
x = 1 - x - 2 * y;
cout << " $$ ( " << x << ", " << y << " ) $$ " << endl;
}
Note that this , well, path, moves in on one spiral arm then back out on the other; they really make up a single path.
$$ ( 1189, 2871 ) $$
$$ ( 1189, -492 ) $$
$$ ( -204, -492 ) $$
$$ ( -204, 85 ) $$
$$ ( 35, 85 ) $$
$$ ( 35, -14 ) $$
$$ ( -6, -14 ) $$
$$ ( -6, 3 ) $$
$$ ( 1, 3 ) $$
$$ ( 1, 0 ) $$
$$ ( 0, 0 ) $$
$$ ( 0, 1 ) $$
$$ ( -1, 1 ) $$
$$ ( -1, -2 ) $$
$$ ( 6, -2 ) $$
$$ ( 6, 15 ) $$
$$ ( -35, 15 ) $$
$$ ( -35, -84 ) $$
$$ ( 204, -84 ) $$
$$ ( 204, 493 ) $$
$$ ( -1189, 493 ) $$
$$ ( -1189, -2870 ) $$
$$ ( 6930, -2870 ) $$
$$ ( 6930, 16731 ) $$
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Moving outward on the black spiral, which contains (1,3)
mpz_class x,y;
x = 1; y = 0;
cout << " $$ ( " << x << ", " << y << " ) $$ " << endl;
while ( abs(x) < 2000000000 && abs(y) < 2000000000)
{
y = 1 + 2 * x - y;
cout << " $$ ( " << x << ", " << y << " ) $$ " << endl;
x = 1 - x - 2 * y;
cout << " $$ ( " << x << ", " << y << " ) $$ " << endl;
}
$$ ( 1, 0 ) $$
$$ ( 1, 3 ) $$
$$ ( -6, 3 ) $$
$$ ( -6, -14 ) $$
$$ ( 35, -14 ) $$
$$ ( 35, 85 ) $$
$$ ( -204, 85 ) $$
$$ ( -204, -492 ) $$
$$ ( 1189, -492 ) $$
$$ ( 1189, 2871 ) $$
$$ ( -6930, 2871 ) $$
$$ ( -6930, -16730 ) $$
$$ ( 40391, -16730 ) $$
$$ ( 40391, 97513 ) $$
$$ ( -235416, 97513 ) $$
$$ ( -235416, -568344 ) $$
$$ ( 1372105, -568344 ) $$
$$ ( 1372105, 3312555 ) $$
$$ ( -7997214, 3312555 ) $$
$$ ( -7997214, -19306982 ) $$
$$ ( 46611179, -19306982 ) $$
$$ ( 46611179, 112529341 ) $$
$$ ( -271669860, 112529341 ) $$
$$ ( -271669860, -655869060 ) $$
$$ ( 1583407981, -655869060 ) $$
$$ ( 1583407981, 3822685023 ) $$
$$ ( -9228778026, 3822685023 ) $$
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Moving outward on the purple spiral, which contains (6,15)
mpz_class x,y;
x = 0; y = 0;
cout << " $$ ( " << x << ", " << y << " ) $$ " << endl;
while ( abs(x) < 2000000000 && abs(y) < 2000000000)
{
y = 1 + 2 * x - y;
cout << " $$ ( " << x << ", " << y << " ) $$ " << endl;
x = 1 - x - 2 * y;
cout << " $$ ( " << x << ", " << y << " ) $$ " << endl;
}
$$ ( 0, 0 ) $$
$$ ( 0, 1 ) $$
$$ ( -1, 1 ) $$
$$ ( -1, -2 ) $$
$$ ( 6, -2 ) $$
$$ ( 6, 15 ) $$
$$ ( -35, 15 ) $$
$$ ( -35, -84 ) $$
$$ ( 204, -84 ) $$
$$ ( 204, 493 ) $$
$$ ( -1189, 493 ) $$
$$ ( -1189, -2870 ) $$
$$ ( 6930, -2870 ) $$
$$ ( 6930, 16731 ) $$
$$ ( -40391, 16731 ) $$
$$ ( -40391, -97512 ) $$
$$ ( 235416, -97512 ) $$
$$ ( 235416, 568345 ) $$
$$ ( -1372105, 568345 ) $$
$$ ( -1372105, -3312554 ) $$
$$ ( 7997214, -3312554 ) $$
$$ ( 7997214, 19306983 ) $$
$$ ( -46611179, 19306983 ) $$
$$ ( -46611179, -112529340 ) $$
$$ ( 271669860, -112529340 ) $$
$$ ( 271669860, 655869061 ) $$
$$ ( -1583407981, 655869061 ) $$
$$ ( -1583407981, -3822685022 ) $$
$$ ( 9228778026, -3822685022 ) $$
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
Best Answer
Let $Y=2t-1$ and $X=2r-1$.
Then $Y^2-2X^2=-1$.
This is the negative Pell equation $y^2-nx^2=-1$ with $n=2$.
Solutions are $Y=$$1, 7, 41, 239,...$
and $X=$$1,5,29,169,...$,
so $t=0,4,21,120,...$
and $r=0,3,15,85,....$
[Click on the numbers to see more solutions from The On-Line Encyclopedia of Integer Sequences.]