Suppose we wish to solve $$x^2-D{y}^2=1$$ where $D$ is not a square. First we compute the continued fraction expansion of $\sqrt{D}$. $$\begin{array}{ccc}{P}_0& =& 0\\ Q_0& =& 1\\ a_0& =& \lfloor (\sqrt{D})\rfloor \\ P_1& =& a_0\\ Q_1& =& D-a_0^2\\ a_n& =& \lfloor \frac{a_0+P_n}{Q_n}\rfloor \\ P_n& =& a_{n-1}{Q}_{n-1}-P_{n-1}\\ Q_n& =& \frac{D-P_n^2}{Q_{n-1}}\end{array}$$

For some $k$ we have $a_{k+1}=2a_0$. Then the continued fraction expansion of $\sqrt{D}$ is $$\sqrt{D}=[a_0{a}_1...a_{k+1}{a}_1...a_{k+1}...]$$

The convergents are computed by $$\begin{array}{ccc}{p}_0& =& a_0\\ p_1& =& a_0{a}_1+1\\ p_n& =& a_n{p}_{n-1}+p_{n-2}\\ q_0& =& 1\\ q_1& =& a_1\\ q_n& =& a_n{q}_{n-1}+q_{n-2}\\ \end{array}$$ The convergents satisfy $$p_n^2-D{q}_n^2=(-1{)}^{n+1}{Q}_{n+1}$$ It turns out $(x,y)=(p_k,q_k)$ is the smallest integer solution of the Pell equation for odd $k$, and $(x,y)=(p_{2k+1},q_{2k+1})$ if $k$ is even. From a minimal positive solution $(t,u)$, we may generate the other positive solutions $(x,y)$ via $$x+y\sqrt{D}=(t+u\sqrt{D}{)}^n$$ for all positive $n$.

Generalized Pell Equations

Now consider the equation $$x^2-D{y}^2=N$$ If $N^2<D$ we may do the following. Compute the convergents $p_n,q_n$ until the smallest integer solution of the Pell equation $x^2-D{y}^2=1$ is found. In the meantime, check if each $p_n,q_n$ satisfies $p_n^2-D{q}_n^2=N/f^2$ for some $f>0$. If so, append $(f{p}_n,f{q}_n)$ to the list of solutions.

We now use this list of solutions to generate all other solutions. If $(r,s)$ is on the list, and $(t,u)$ is a minimal positive solution of the corresponding Pell equation, then we have a family of solutions $(x,y)$ given by $$x+y\sqrt{D}=(r+s\sqrt{D})(t+u\sqrt{D}{)}^n$$ for all positive $n$.

If $N^2\ge D$ then one possibility is brute force. If $N>0$ set $L_1=0,L_2=\sqrt{N(t-1)/2D)}$, otherwise set $L_1=\sqrt{-N/D},L_2=\sqrt{-N(t+1)/(2D)}$. For $L_1\le y\le L_2$ check if $N+D{y}^2=x^2$ for some integer $x$. If so, then append $(x,y)$ to the list of fundamental solutions. Also append $(-x,y)$ if not equivalent to $(x,y)$. (TODO: can avoid since we don't want negative solutions?) Then proceed as above to generate all solutions.