Let $u$ be an element of $\mathbb{Z}[\sqrt 5]$ of norm 1, i.e. $u = r + s \sqrt 5$ with $r^2-5s^2 = 1$.
The multiplication by $u$ in $\mathbb{Z}[\sqrt 5]$ turns any element $y$ of norm $44$ into another element $uy$ of norm $44$.
View this multiplication operation on $\mathbb{Z}[\sqrt 5]$ as the transformation of the plane $f : (p,q) \rightarrow (pr+5qs,ps+qr)$, and look for its eigenvalues :
$f(\sqrt5,1) = (r\sqrt5+5s,r+\sqrt5s) = u(\sqrt5,1)$, and we have $f(- \sqrt5,1) = \frac 1u (- \sqrt5,1)$ as well.
If $u>1$ this means that $f$ is an operation that, when iterated, takes elements near the line $(p = - \sqrt5 q)$ and moves them over to the line $(p = \sqrt5 q)$
Now you want to find a sector of the plane so that you can reach the whole plane by taking its images by the iterates of $f$ and $f^{-1}$
Define $g(p,q) = \frac {p + \sqrt5 q}{p - \sqrt5 q}$, which is the ratio of the coordinates of $(p,q)$ in the eigenbasis of $f$.
$g(f(p,q)) = \frac {pr+5qs + \sqrt5 (ps+qr)}{pr+5qs - \sqrt5 (ps+qr)} = \frac{(r+\sqrt5 s)(p + \sqrt5 q)}{(r-\sqrt5 s)(p - \sqrt5 q)} = (r+\sqrt5 s)^2 g(p,q)$.
Or alternately, define $g(y) = y/\overline{y}$, so that $g(uy) = uy/\overline{uy} = u^2 g(y)$.
Thus if you look at any point $(p,q)$, you know you can apply $f$ or $f^{-1}$ to turn it into $(p',q')$ such that $g(p',q') \in [1 ; u^2[$
Thus, a suitable sector of the plane is the set of points $(p,q)$ such that $g(p,q) \in [1 ; u^2[$ : if you find all the elements $y$ of norm $44$ such that $g(y) \in [1 ; u^2[$, then this means that the $u^ky$ will cover all the elements of norm $44$
Finally, the good thing is that $ \{y \in \mathbb{Z}[ \sqrt 5] / g(y) \in [1; u^2[, y\overline{y} \in [0; M] \}$ is a finite set, so a finite computation can give you all the elements of norm $44$ you need.
In the case of $p²-10q²=9$, a fundamental unit is $u = 19+6\sqrt{10}$,
so replace $\sqrt 5,r,s$ with $ \sqrt {10},19,6$ in everything I wrote above.
In order to find all the solutions, you only need to check potential solutions in the sector of the plane between the lines $g(p,q) = 1$ and $g(p,q) = u^2$.
You can look at the intersection of the line $g(p,q)=1$ with the curve $p^2-10q^2 = 9$.
$g(p,q)=1$ implies that $p+\sqrt{10}q = p- \sqrt{10}q$, so $q=0$, and then the second equations has two solutions $p=3$ and $p= -3$.
It so happens that the intersection points have integer coordinates so they give solutions to the original equation.
Next, the intersection of the line $g(p,q) = u^2$ with the curve will be $u \times (3,0) = f(3,0) = (19*3+60*0, 6*3+19*0) = (57,18)$ and $u \times (-3,0) = (-57,-18)$.
So you only have to look for points on the curve $p^2-10q^2=9$ with integers coordinates in the section of the curve between $(3,0)$ and $(57,18)$ (and the one between $(-3,0)$ and $(-57,-18)$ but it is essentially the same thing).
You can write a naïve program :
for q = 0 to 17 do :
let square_of_p = 9+10*q*q.
if square_of_p is a square, then add (sqrt(square_of_p),q) to the list of solutions.
Which will give you the list $\{(3,0) ; (7,2) ; (13,4)\}$.
These three solutions, together with their opposite, will generate, using the forward and backward interations of the function $f$, all the solution in $\mathbb{Z}^2$.
If you only want solution with positive coordinates, the forward iteration of $f$ on those three solutions are enough.
Also, as Gerry points out, the conjugate of $(7,2)$ generates $(13,4)$ because $f(7,-2) = (13,4)$. Had we picked a sector of the plane symmetric around the $x$-axis, we could have halved the search space thanks to that symmetry, and we would have obtained $\{(7,-2),(3,0),(7,2)\}$ instead.
One loop of this hypnotic animation represents one application of the function $f$.
Each dot corresponds to one point of the plane with integer coordinates, and is moved to its image by $f$ in the course of the loop.
The points are colored according to their norm (and as you can see, each of them stay on their hyperbolic branch of points sharing their norm), and I've made the yellow-ish points of norm 9 (the solutions of $x^2-10y^2 = 9$) a bit bigger.
For example, the point at (3,0) is sent outside the graph, and the point at (-7,2) is sent on (13,4) (almost vanishing).
You can see that there are three points going through (3,0) during the course of one loop. They correspond to three representants of the three fundamental solutions of the equation.
For each yellowish point on the curve $x^2-10y^2=9$, no matter how far along the asymptote it may be, there is an iterate of $f$ or $f^{-1}$ that sends it to one of those three fundamental solutions.
In order to find all fundamental solutions, it is enough to explore only a fundamental portion of the curve (a portion whose iterates by $f$ covers the curve), for example the fundamental portion of the curve between (-7,2) and its image by $f$, (13,4). To find the solutions on that portion, you set $y=-2,-1,0,1,2,3$ and look if there is an integer $x$ that makes a solution for each of those $y$.
Whichever fundamental portion of the curve you choose, you will find 3 solutions inside it, whose images by $f$ are sent to the next three solutions in the next portion of the curve, and so on.
Now there is a better procedure than the "brute search" I did to get all the solutions.
It is an adaptation of the procedure to obtain a fundamental unit :
Start with the equation $x^2-10y^2 = 9$, and suppose we want all the positive solutions.
We observe that we must have $x > 3y$, or else $-y^2 \ge 9$, which is clearly impossible.
So, replace $x$ with $x_1 + 3y$.
We get the equation $x_1^2 + 6x_1 y - y^2 = 9$.
We observe that we must have $y > 6x_1$, or else $x_1^2 \le 9$.
In this case we quickly get the three small solutions $(x_1,y) = (1,2),(1,4),(3,0)$ which correspond to the solutions $(x,y) = (7,2),(13,4),(3,0)$.
Otherwise, continue and replace $y$ with $y_1 + 6x_1$.
We get the equation $x_1^2 - 6x_1y_1 - y_1^2 = 9$.
We observe that we must have $x_1 > 6y_1$, or else $-y_1^2 \ge 9$, which is clearly impossible.
So, replace $x_1$ with $x_2 + 6y_1$.
We get the equation $x_2^2 + 6x_2y_1 - y_1^2 = 9$.
But we already encountered that equation so we know how to solve it.
Best Answer
$$ 8x^2 - 4 y^2 - 4 y - 1 = -1 $$ $$ (2y+1)^2 - 8 x^2 = 1 $$ $$ w^2 - 8 x^2 = 1 \; , $$ so that $w$ really is odd, then $y = (w-1)/2$
The first two are $$ (w,x) = (1,0) $$ then $$ (w,x) = (3,1) $$
After that, growth comes from $$ (w,x) \mapsto (3w+8x, w + 3x ) $$ over and over
but still followed by $y = (w-1)/2$ for each
If preferred, there are separate recurrences
$$ w_{n+2} = 6 w_{n+1} - w_n \; , \; $$ $$ x_{n+2} = 6 x_{n+1} - x_n \; . \; $$