Hello Akiva,
The Diophantine equation you are considering is closely related to
Pell's equation.
You have made the astute observation that the solutions to this equation are related to Pythagorean triples where the legs differ by 1. You may find http://www.mathhelpboards.com/f49/pell-sequence-2905/ to be of interest.
Okay, you have found the solutions:
$(x,y)=(1,1),\,(7,5),\,(41,29)$
At this point, we could hypothesize that $x$ and $y$ can be given recursively with:
$$S_{n+1}=6S_{n}-S_{n-1}$$
which leads to the characteristic equation:
$$r^2-6r+1=0$$
which has the roots:
$$r=3\pm2\sqrt{2}=\left(1\pm\sqrt{2} \right)^2$$
and so we know the closed form is:
$$S_n=k_1\left(1-\sqrt{2} \right)^{2n}+k_2\left(1+\sqrt{2} \right)^{2n}$$
Next, we may using the initial values for each sequence to determine the parameters.
For $x_n$, we may use:
$$x_1=k_1\left(1-\sqrt{2} \right)^{2}+k_2\left(1+\sqrt{2} \right)^{2}=1$$
$$x_2=k_1\left(1-\sqrt{2} \right)^{4}+k_2\left(1+\sqrt{2} \right)^{4}=7$$
This gives rise to the system:
$$3\left(k_1+k_2 \right)+2\sqrt{2}\left(k_2-k_1 \right)=1$$
$$17\left(k_1+k_2 \right)+12\sqrt{2}\left(k_2-k_1 \right)=7$$
Multiplying the first equation by $-6$, then adding, we get:
$$-k_1-k_2=1\,\therefore\,k_1=-1-k_2$$
Substituting into the first equation for $k_1$, there results:
$$-3+2\sqrt{2}\left(2k_2+1 \right)=1$$
Solving for $k_2$, we find:
$$k_2=\frac{\sqrt{2}-1}{2}$$
Thus:
$$k_1=-1-\frac{\sqrt{2}-1}{2}=-\frac{\sqrt{2}+1}{2}$$
and so we find:
$$x_n=-\frac{\sqrt{2}+1}{2}\left(1-\sqrt{2} \right)^{2n}+\frac{\sqrt{2}-1}{2}\left(1+\sqrt{2} \right)^{2n}$$
$$x_n=\frac{1}{2}\left(\left(\sqrt{2}-1 \right)\left(1+\sqrt{2} \right)^{2n}-\left(\sqrt{2}+1 \right)\left(1-\sqrt{2} \right)^{2n} \right)$$
For $y_n$, we may use:
$$y_1=k_1\left(1-\sqrt{2} \right)^{2}+k_2\left(1+\sqrt{2} \right)^{2}=1$$
$$y_2=k_1\left(1-\sqrt{2} \right)^{4}+k_2\left(1+\sqrt{2} \right)^{4}=5$$
This gives rise to the system:
$$3\left(k_1+k_2 \right)+2\sqrt{2}\left(k_2-k_1 \right)=1$$
$$17\left(k_1+k_2 \right)+12\sqrt{2}\left(k_2-k_1 \right)=5$$
Multiplying the first equation by $-6$, then adding, we get:
$$-k_1-k_2=-1\,\therefore\,k_1=1-k_2$$
Substituting into the first equation for $k_1$, there results:
$$3+2\sqrt{2}\left(2k_2-1 \right)=1$$
Solving for $k_2$, we find:
$$k_2=\frac{\sqrt{2}-1}{2\sqrt{2}}$$
Thus:
$$k_1=1-\frac{\sqrt{2}-1}{2\sqrt{2}}=\frac{\sqrt{2}+1}{2\sqrt{2}}$$
and so we find:
$$y_n=\frac{\sqrt{2}+1}{2\sqrt{2}}\left(1-\sqrt{2} \right)^{2n}+\frac{\sqrt{2}-1}{2\sqrt{2}}\left(1+\sqrt{2} \right)^{2n}$$
$$y_n=\frac{1}{2\sqrt{2}}\left(\left(\sqrt{2}+1 \right)\left(1-\sqrt{2} \right)^{2n}+\left(\sqrt{2}-1 \right)\left(1+\sqrt{2} \right)^{2n} \right)$$
So, in summary, we find an infinite number of solution pairs, given:
i) recursively:
$$x_{n+1}=6x_{n}-x_{n-1}$$ where $$x_1=1,\,x_2=7$$
$$y_{n+1}=6y_{n}-y_{n-1}$$ where $$y_1=1,\,x_2=5$$
ii) closed form:
$$x_n=\frac{1}{2}\left(\left(\sqrt{2}-1 \right)\left(1+\sqrt{2} \right)^{2n}-\left(\sqrt{2}+1 \right)\left(1-\sqrt{2} \right)^{2n} \right)$$
$$y_n=\frac{1}{2\sqrt{2}}\left(\left(\sqrt{2}+1 \right)\left(1-\sqrt{2} \right)^{2n}+\left(\sqrt{2}-1 \right)\left(1+\sqrt{2} \right)^{2n} \right)$$