Solution 1:

$$(\frac{x}y)^2-(\frac{1}y)^2=D,\\(\frac{x}y+\frac{1}y)(\frac{x}y-\frac{1}y)=D,\\\frac{x}y+\frac{1}y=t,\frac{x}y-\frac{1}y=\frac{D}t,\\\frac{x}y=\frac{1}2(t+\frac{D}{t}),\frac{1}y=\frac{1}2(t-\frac{D}{t}),\\x=\dfrac{t^2+D}{t^2-D},y=\frac{2t}{t^2-D},(t\in \mathbb Q,t^2\neq D).$$

Solution 2:

The following theorem is due to David Hilbert:

Let $K/\mathbf Q$ be a cyclic Galois extension, with Galois group generated by $\sigma$. Then every element of norm $1$ of $K$ is of the form $\sigma(x)/x$ for some $x \in K^\times$.

Applying this to $K=\mathbf Q(\sqrt D)$ and noticing that $N(x+\sqrt D y) = x^2-Dy^2$, it follows that every solution to $N(x+\sqrt D y) = 1$ is of the form

$$x+\sqrt Dy = \frac{u - \sqrt D v}{u + \sqrt D v} = \frac{u^2+Dv^2}{u^2-Dv^2} + \sqrt D \frac{-2uv}{u^2-Dv^2},$$

i.e.

$$(x,y) = \left(\frac{u^2+Dv^2}{u^2-Dv^2}, \frac{-2uv}{u^2-Dv^2}\right), \qquad (u, v) \in \mathbf Q,\qquad u^2+v^2 \neq 0.$$

Solution 3:

There is a simple and more general way to look at this problem that applies to finding rational points on any conic section with rational coefficients. i.e. any equation in the form $ax^2 + bxy + cy^2 + dx + ey + f = 0$ , $a,b,c,d,e,f \in Q$

The Pell equation is just one example.

The solution is to find one rational point and then draw a straight line through that point at some rational valued gradient $m$. Wherever this line meets the conic will be another rational point and all rational points can be found this way from just the one fixed rational point.

The reason this works is that when you substitute the equation for the line into the conic you get a quadratic equation, so if you have one rational root the other root of the quadratic will be another rational root.

Let's see how this works for Pell's equation. $x^2 - Dy^2 = 1$

One rational root is given by $y=0, x=1$. The line through this point with gradient $m$ is $y=m(x-1)$. Substitute into Pell's equation to find the second intersect.

$x^2 - D(m(x-1))^2 = 1$

$(x-1)(x+1-Dm^2(x-1)) = 0$

So the second rational point is at

$x+1-Dm^2(x-1) = 0$

$Dm^2+1 = x(Dm^2-1)$

$x = \frac{Dm^2+1}{Dm^2-1}$

Use the equation for the line to get the corresponding value of $y$

$y=m(x-1) = \frac{2m}{Dm^2-1}$

Note that you should also use the value $m$ at infinity which in this case gives back the point $(x,y)=(1,0)$