Let Peel's Equation be
$$x^2-D\cdot y^2=1.\tag{P}$$
where $D$ is a given positive nonsquare integer.
In solving (P) we must find the fundamental solution, i.e. a pair of positive integers $(x_1,y_1)$ wich satisfie them,and $x_1$ is minimal. An algorithm that determines this, found by me in the book "Dr. Paul RADOVICI-MĂRCULESCU, Probleme de Teoria Elementară a Numerelor" (Ed. Tehnică, București, 1986), at page 178, uses Continued Fraction.
We develop $\sqrt{D}$ in continued fraction (their shape is well known):
$\sqrt{D}=[a_0,\; \overline{a_1,\;\dots,\;a_{k-1},\;2a_0}] \tag{D}$
$\blacklozenge$ if $k$=odd, $k=2t+1$ then the fundamental solution is $x_1=A_{4t+1},\;y_1=B_{4t+1}; \tag{O}$
$\blacklozenge$ if $k$=even, $k=2t$ then the fundamental solution is $x_1=A_{2t-1},\;y_1=B_{2t-1}.\tag{E}$
I wrote $A_k$ the numerator and $B_k$ the denominator of the $k$th convergent of (D), i.e.
$[a_0,\;a_1,\;\dots, a_k]=\frac{A_k}{B_k}.$
I found their calculation formulas in the book ""I.M. VINOGRADOV, Bazele Teoriei Numerelor"(Ed. Academiei, 1957), at page 17
$$A_i=a_i\cdot A_{i-1}+A_{i-2},\;\;B_i=a_i\cdot B_{i-1}+B_{i-2},\;\;i\geqslant 1$$
with initial values $A_{-1}=1,\;B_{-1}=0;\,\,\,A_0=a_0,\;B_0=1$
where An is the numerator and Bn is the denominator, called continuants, of the nth convergent.
Schematically, obtaining $A_i$ can be represented as follows
$\begin{cases}\;\;\;\;\;\;\;\;\times \swarrow a_i\\\underset{+}{\underbrace{A_{i-2}\;A_{i-1}}}\;\fbox{Ai}\tag{S}\end{cases}$
And similar for $B_i$.
Example 1 $x^2-13\cdot y^2=1$
We have the continued fraction expansion
$\sqrt{13}=[3,\;\overline{1,1,1,1,6}]$ so $k=5$.
We apply the scheme (S), the calculated elements being written in color, and we obtain the table:
\begin{array}{|c|c}\hline\\\;i&-1&0&1&2&3&4&5&6&7&8&\color{Green}9&10\\\hline \\a_i&\;&3&1&1&1&1&6&1&1&1&1&6\\\hline \\A_i&1&3&\color{Red}4&\color{Red}7&\color{Red}{11}&\color{Red}{18}&\color{Red}{119}&\color{Red}{137}&\color{Red}{256}&\color{Red}{393}&\color{Green}{649}&\color{Red}{4287}\\\hline\\B_i&0&1&\color{Red}1&\color{Red}2&\color{Red}3&\color{Red}5&\color{Red}{33}&\color{Red}{38}&\color{Red}{71}&\color{Red}{109}&\color{Green}{180}&\color{Red}{1189}\\\hline \end{array}
We are in the case (O), $k=5=2\cdot t+1\;\Rightarrow t=2$ and the minimal solution is (colored green)
$x_1=A_9=649,\;\;y_1=B_9=180.$
Remark CiP The result is in agreement with that given in the Table on page 127 of the book
ANDREESCU Titu, ANDRICA Dorin, CUCUREZEANU Ion - An Introduction to Diophantine Equations : A Problem-Based Approach (Birkhäuser Springer, New York Dordrecht Heidelberg London, 2010)
<end Rem>
Example 2 $x^2-31\cdot y^2=1$
We have the continued fraction expansion
$\sqrt{31}=[5,\;\overline{1,1,3,5,3,1,1,10}]$ so $k=8.$
Now we are in the case (E): $8=k=2\cdot t\;\Rightarrow\;t=4$ and the minimal solution is, according to the calculations in the table below
$x_1=A_7=1520,\;\;y_1=B_7=273.$
(Compare with the Table mentioned in the previous Remark.)
\begin{array}{|c|c|}\hline\\i&-1&0&1&2&3&4&5&6&\color{Green}7\\\hline \\a_i&\;&5&1&1&3&5&3&1&\color{Green}1\\\hline\\A_i&1&5&6&11&39&206&657&863&\color{Green}{1520}\\\hline\\B_i&0&1&1&2&7&37&118&155&\color{Green}{273}\\\hline\end{array}
As for ALL the solutions in natural numbers of the equation (P), they are given by the formulas (O) - wiyh $t=0,\,1,\,2,\dots$ and (E) - with $t=1,\;2,\dots$. However, the further calculation of $n$th convergent is too laborious.
A simpler approach involves writing a system of recurrence relations.
Suppose, for the equation (P) that we have found the fundamental solution $(x_1,y_1)$. Then, using the matrix
$P=P_D=\begin{pmatrix}x_1&D\cdot y_1\\y_1&x_1\\\end{pmatrix} \tag{M}$
we write the linear recurrence relation in the form
$\begin{pmatrix}x_{n+1}\\y_{n+1} \end{pmatrix}=P \cdot \begin{pmatrix}x_n\\y_n \end{pmatrix}=\begin{pmatrix}x_1&D\cdot y_1\\y_1&x_1 \end{pmatrix} \cdot \begin{pmatrix}x_n\\y_n \end{pmatrix}\;,n=0,\;1,\dots\tag{MX}$
with $(x_0,y_0)=(1,0)$ which is the universal solution for the equations (P),
that is
$\begin{cases}x_{n+1}=x_1\cdot x_n+(Dy_1)\cdot y_n\\y_{n+1}=y_1\cdot x_n+x_1 \cdot y_n \end{cases},\;\;x_0=1,\;y_0=0,\;\;\;n=0,\;1,\dots \tag{XY}$
The two-recurrence system (XY) can be transformed into second-order recurrence relations, as shown in another post.
Example 3 $x^2-3\cdot y^2=1$
Fundamental solution of this equation is $(x_1,y_1)=(2,1)$. So $P=\begin{pmatrix}2&3\cdot 1\\1&2\end{pmatrix}$ and the recurrence relations are written
$\begin{cases}x_{n+1}=2\cdot x_n+3\cdot y_n\\y_{n+1}=x_n+2 \cdot y_n\end{cases},\;x_0=1,\;y_0=0\;\;\;(n=0,\;1,\dots).$
The first few solutions appear in the Table below:
\begin{array}{c|c}n&0&1&2&3&4&5&6\\\hline \\x_n&1&2&7&26&97&362&1351\\\hline \\y_n&0&1&4&15&56&209&780 \end{array}
Or we can write the second-order recurrence relations, as mentioned above, with the same results:
$x_{n+2}=4\cdot x_{n+1}-x_n,\;\;x_0=1,\;x_1=2\;\;\;\;;\;\;\;\;y_{n+2}=4\cdot y_{n+1}-y_n,\;\;y_0=0,\;y_1=1\;\;(n=0,\;1,\dots)$.
Now let the Negative Pell (or "minus Pell") equation be
$x^2-D\cdot y^2=-1 \tag{-P}$
Not all of these equations have solutions. But when they do, they solve similarly. I first encountered it in an article in GM-P, no. 4/1989, pages 175-178: Gheorghe UDREA - "Asupra Formei Solutiilor Ecuatiilor $x^2-D\cdot y^2=\pm1$". We quote from here what interests us.
Let's go back to (D).
$\quad\blacklozenge$ if $k$=odd the solutions of the equation (-P) are
$x_n=A_{k(2n-1)-1},\;\;\;y_n=B_{k(2n-1)-1},\;\;\;n=1,\;2,\;3,\dots ;\tag{-O}$
$\quad\blacklozenge$ if $k$=even the equation (-P) has NO solutions.
So, the smallest positive solution, when it exists (so $k$ is odd) is
$x_1=A_{k-1}\,,\,\,\,y_1=B_{k-1} \tag{-F}$
Furthermore, we can also determine the smallest positive solution(the fundamental one) of the equation (P) by the formula
$x_1'+y_1'=(x_1+y_1\cdot \sqrt{D})^2 \tag{+P}$
Example 4 $x^2-13\cdot y^2=-1$
We have (see Example 1) $\sqrt{13}=[3,\overline{1,1,1,1,6}]$ so $k=5$. Minimal solution is $x_1=A_4=18,\;y_1=B_4=5.$ Further we have
$(18+5\cdot \sqrt{13})^2=324+25\cdot 13+180\cdot \sqrt{13}=649+180 \cdot \sqrt{13}$
so minimai solution of $x^2-13\cdot y^2=1$ is $x_1'=649,\;y_1'=180$ as we have already found.
Images of the older manuscript and some drafts.