We'll study here only the case D positive not a square and first of all
the case k = 1 :
x²  Dy² = 1 (proper Pell equation).
We then have the following results :
x_{n+1} = x_{0}x_{n} + Dy_{0}y_{n}, y_{n+1} = y_{0}x_{n} + x_{0}y_{n}
We can then deduce recurence relation for the sequences x_{n} and y_{n} :
x_{n+2} = 2x_{0}x_{n+1}x_{n} y_{n+2} = 2x_{0}y_{n+1}y_{n}
We can also write : x_{n} + y_{n}√D = (x_{0} + y_{0}√ D)^{n+1} or the conjugate form, replacing + by 
This gives closed forms for x_{n} and y_{n} as a function of n.
These formulaes seems to be of little use because of irrational terms in √D, although giving
integral x_{n} and y_{n} !
(They are easy to write from the previous relation combined with its conjugate)
Let U_{0} = 0, V_{0} = 1, a_{0} = E[√D], P_{0} = a_{0}, Q_{0} = 1
and by convention P_{1} = 1 and Q_{1} = 0
Noting E[x] = integer part of x. Then repeat : U_{n+1} = a_{n}V_{n}  U_{n},
V_{n+1} = (D  U_{n+1}²)/V_{n},
a_{n+1} = E[ (a_{0} + U_{n+1})/V_{n+1} ],
until a_{k} > a_{0} (or V_{k} = 1).
If k even it is the required solution.

This algorithm is a specific case of PQa algorithm giving the expansion as a continued fraction of a quadratic number (u + √D)/v with here u = 0 and v = 1.
Apply this to x²  13y² = 1 :
U_{0} = 0  V_{0} = 1  a_{0} = 3  P_{0} = 3  Q_{0} = 1  
U_{1} = 3  V_{1} = (13  9)/1 = 4  a_{1} = E[(3 + 3)/4] = 1  P_{1} = 4  Q_{1} = 1  
U_{2} = 4  3 = 1  V_{2} = (13  1)/4 = 3  a_{2} = E[(3 + 1)/3] = 1  P_{2} = 7  Q_{2} = 2  
U_{3} = 3  1 = 2  V_{3} = (13  4)/3 = 3  a_{3} = E[(3 + 2)/3] = 1  P_{3} = 11  Q_{3} = 3  
U_{4} = 3  2 = 1  V_{4} = (13  1)/1 = 4  a_{4} = E[(3 + 1)/4] = 1  P_{4} = 18  Q_{4} = 5  
U_{5} = 4  1 = 3  V_{5} = (13  9)/4 = 1  a_{5} = E[(3 + 3)/1] = 6 > 3 => stop 
5 being odd, (18, 5) is solution of x²  13y² = 1.
The fundamental solution of x²  13y² = 1 is (18² + 13×5², 2×18×5) = (649, 180)
which was not obvious without this algorithm.
Let's go on with the general case k ≠ 1, and first with
or the conjugate form by replacing + by 
Note : solutions of equation x²  Dy² = 1 are obtained from the fundamental solution
(x_{0}, y_{0}) of x²  Dy² = 1
by : x_{n} + y_{n}√D = (x_{0} + y_{0}√D)^{2n}
The fundamental solution of x²  Dy²=1 is then (x_{0}² + Dy_{0}², 2x_{0}y_{0}).
We also have the recurrence relations :
x_{n+1} = (x_{0}² + Dy_{0}²)x_{n} + 2Dx_{0}y_{0}y_{n}, y_{n+1} = 2x_{0}y_{0}x_{n} + (x_{0}² + Dy_{0}²)y_{n}
As well as relations in x_{n+2} = ax_{n+1} + bx_{n} that are left as an exercice.
As for the searching of the fundamental solution, just use the previous algorithm, and note that when k is even,
then x²  Dy² = 1 has no solution.
There is no easy criterion to know it in advance. Note the case when D is a prime p :
Equation x²  py² = 1 has solutions if and anly if p = 4k + 1.
Equation x²  2y² = 1
x = 1, y = 1 is an obvious solution, and is fundamental as any smaller solutions would have
x and y < 1.
The fundamental solution of x²  2y² = +1 is then x = 1² + 2*1² = 3
and y = 1² + 1² = 2 .
All other solutions are then
x_{n+1} = 3x_{n} + 4y_{n}, y_{n+1} = 2x_{n} + 3y_{n}, that is :
(1,1) (7,5) (41,29) (239,169) (1393,985) ...
x_{n} + y_{n}√D = (x_{0} + y_{0}√D)^{n}(ξ + η√D) and also the corresponding recurrence relations.
The problem is now to find these independant primitive solutions. A first method is to scan systematically all values of x or y, the theory giving just an upper bound. (remind the equation may have no solutions). This is the most efficient method if the upper bound is not too big.
But at first we must explain what are independant solutions,
or more precisely when they are not, that is when they are said to be "equivallent".
Two solutions (x,y) and (x',y') are equivalent if there is a solution (X,Y) of
x²  Dy² = 1 with :
x + y√D = (X + Y√D)(x' + y'√D)
This criterion being not very usefull, add the following :
Two solutions (x,y) and (x',y') of x²  Dy² = K are equivallent
if and only if :
xx'  Dyy' and x'y  xy' are multiple of K. 
Remains to find bounds for the primitive solutions :
If K > 0, the primitive solutions are bounded by
0 < x ≤ √( (x_{0} + 1)K/2 ) or 0 ≤ y ≤ y_{0}√ ( K/2(x_{0} + 1) ) If K < 0, the primitive solutions are bounded by

If there is no solutions in the bounds, there is no solutions at all.
We would like to find a simple criterion for x²  Dy² = K to have solutions by writing
x² ≡ K mod (D),
that is K must be a quadratic residue of D.
The only thing we can say then is that when K is not a quadratic residue of D,
then the Pell equation has no solutions.
For instance equation x²  97y² = 51 has no solutions Legendre symbol (51/97) = 1
hence 51 is not a quadratic residue of 97.
But this necessary condition is not sufficient...
Hence there are 4 primitive solutions, but not all are independant
(±3, 0) are equivalent, more generally (x,y) ~ (x,y)
(4,±1) are independant : 16 + 7 = 26 coprime with 9.
(3,0) and (4,1) are independant
(3,0) and (4,1) are independant
Finally three independant primitive solutions : (3, 0), (4, 1) and (4,1)
If we just care of x,y > 0, (4,1) ~ (11,4), is then smallest >0 solution in this family
All solutions are built from the three families :
x_{n} + y_{n}√7 = (8 + 3√7)^{n}(3 + 0√7)
x_{n} + y_{n}√7 = (8 + 3√7)^{n}(4 + √7)
x_{n} + y_{n}√7 = (8 + 3√7)^{n}(4  √7)
Or the recurrence relations :
x_{n+1} = 8x_{n} + 21y_{n}, y_{n+1} = 3x_{n} + 8y_{n}, from (x, y) = (3,0), (4,1) and (4,1)
that is the three families :
(3,0) (24,9) (381,144) ...
(11,4) (172,65) (2741,1036) ...
(4,1) (53,20) (844,319) ...
Scanning for solutions becomes tedious if the bound is too big. In that case we have two efficient methods.
P_{i}²  DQ_{i}² = (1)^{i+1}V_{i+1}V_{0}
This allows to find, with the same algorithm, the primitive solutions of x²  Dy² = K
First of all note that if (x,y) is solution of x²  Dy² = K, (mx,my) is a solution of x²  Dy² = m²K.
If we get a (1)^{i+1}V_{i+1} = K/m²,
then x = mP_{i}, y = mQ_{i} is solution of x²  Dy² = K
If K² < D, this gives all primitive solutions.
Example x²  19y² = 3
U  V  a  P  Q  
0  1  4  4  1  1  
4  3  2  9  2  3  4  1  is solution of x²  19y² = 3 
2  5  1  13  3  5  
3  2  3  48  11  2  
3  5  1  61  14  5  
2  3  2  170  39  3  61  14  is solution of x²  19y² = 3 
4  1  8  1  170  39  is solution of x²  19y² = 1 
Notes :
1) Equation x²  19y² = K hence has no solutions for K = 1, +2, +3, and 4.
2) x²  19y² = 5 has primitive solutions (9,2) and (48,11) but we don't know if there are others because 5² > 19.
List the f with f² dividing K, m=K/f²
For each such f, find all z : m/2 < z≤m/2 and z² ≡ D modulo m For each z, apply the PQA algorithm with U_{0} = z, V_{0} = m up to end of period or a V_{i} = ±1  If V_{i} = 1, then (P_{i1},Q_{i1}) is solution of x²  Dy² = m, add this solution to the list.  If V_{i} = 1, then (P_{i1},Q_{i1}) is solution of x²  Dy² = m. If x²  Dy² = 1 has a solution (r,s) (primitive) then (P_{i1},Q_{i1})o(r,s) is solution of x²  Dy²=m. Add (P_{i1},Q_{i1})o(r,s) = (P_{i1}r + Q_{i1}sD, P_{i1}s + Q_{i1}r) 
All these solutions being independant, hence result into 6 solution families by combining them with the solution (649,180)^{n} of x²  13y² = 1