Assuming that $a = x^2$ and $b = y^2$, i converted this equation to a linear diophantine equation for sake of convenience:
$$4a – 3b = 1$$
where after calculating a particular solution (like $(1, 1)$) by extended euclidean algorithm, the $n$th pair of solutions can be calculated as
$$a_n = 1 – 3n$$
$$b_n = 1 – 4n$$
where $n \le 0$ since $a_n, b_n$ are equal to squared numbers and so they have to be equal or greater than zero. So let the:
$$ 0 <= m = -n$$
$$a_m = 3m + 1$$
$$b_m = 4m + 1$$
Then for a pair of $(a_m, b_m)$ where $a_m, b_m$ are perfect squares, the pair $(\sqrt{a_m}, \sqrt{b_m})$ will be a solution for our diophantus equation. Written more generally:
$$(\sqrt{a_m}, \sqrt{b_m}) = (\sqrt{3m + 1}, \sqrt{4m + 1})$$
$$\exists m\in \mathbb Z$$
Of course only when both $3m + 1$ and $4m + 1$ are perfect squares. And this is the point i came at. I don't know what to do after that. Any help would be appreciated.
By the way i wrote a Python script to compute solutions by the method i mentioned above. Here are few pair i succeed to calculate:
$${(1.0, 1.0), (13.0, 15.0), (181.0, 209.0), (2521.0, 2911.0), …}$$
I can't see any connection among the pairs… How i can continue?
Best Answer
As noted, this is the same as $x'-3y^2=1$, known as Pell's equation ($x'=2x$). To find the first nontrivial solution (aka $1,0)$, we take $x'=2,y=1$. Then all solutions are of the form
$$x'=\frac{(2+\sqrt{3})^n+(2-\sqrt{3})^n}{2}$$
$$y=\frac{(2+\sqrt{3})^n-(2-\sqrt{3})^n}{2\sqrt{3}}$$
And $x'$ is even. This occurs precisely when $n$ is odd, as the only terms that are left over in the binomial expansion for the numerator $\mod 4$ are $2n\sqrt{3}^{n-1}$ for $n$ odd, $2\sqrt{3}^n$ for $n$ even.
Thus , we take $x'=2,y=1$. Then all solutions are of the form
$$(x,y)=\left(\frac{(2+\sqrt{3})^{2n-1}+(2-\sqrt{3})^{2n-1}}{2},\frac{(2+\sqrt{3})^{2n-1}-(2-\sqrt{3})^{2n-1}}{2\sqrt{3}}\right)$$
is the final solution set as $n$ ranges over $\mathbb{N}$