Solve in $\mathbf{N}$ the equation $9x^2+p=y^2$ In these days, I have been trying to solve this problem:

Let $p \in \mathbf{N}$ a positive large integer ($> 10^9$). Find all $x, y \in\mathbf{N}$ such that:
$$9x^2+p=y^2$$

The first approach that I have tried is the following. We know that:
$$(t+n)^2-t^2=t^2+2\cdot t \cdot n +n^2-t^2=2\cdot t \cdot n +n^2$$
Also, we can build every possible square $w_n$ greater than $p$ in this way:
$$w_n=\left\lfloor\sqrt{p}\right\rfloor^2+2\cdot \left\lfloor\sqrt{p}\right\rfloor \cdot n +n^2 = \left(n+\left\lfloor\sqrt{p}\right\rfloor^2\right)^2$$
For example, let $p=5$. Follows that: $w_1=\left(1+\left\lfloor\sqrt{5}\right\rfloor^2\right)^2=9$, $w_2=\left(2+\left\lfloor\sqrt{5}\right\rfloor^2\right)^2=16$ and so on.
Now, using this idea and applying to the first equation:
$$9x^2=\left(n+\left\lfloor\sqrt{p}\right\rfloor^2\right)^2-p$$
I am not allowed to apply Pell's equation because $9=3^2$ and calcultaing $\Delta$ in $x$ doesn't help anymore.
Another approach is based on Pell's equation. I thought to express $9x^2=8x^2+x^2$ and then:
$$8x^2+x^2+p=y^2\leftrightarrow 8x^2+p=y^2-x^2\leftrightarrow (y^2-x^2)-8x^2=p \leftrightarrow u^2-8x^2=p$$
But then, in order to generate all the solutions, I have to guess the first one (or one of them) that is pretty complicated for big $p$.
So, how can we do that? Are there any other solutions?
Thanks.
 A: First, reformulate the problem as the following:
$$
y^2 - 9x^2 = p \Rightarrow(y-3x)(y+3x)= p
$$
Now, for any given $p$, find its prime factors. Then, for any 2-partitions of them, solve a simple equation system.
To simplify some cases, suppose $p$ is factorized to $p_1 \times p_2$. Now, solve the following system:
$$
y-3x = p_1
$$
$$
y+3x = p_2
$$
So, $y = \frac{p_1 + p_2}{2}$, and $x = \frac{p_2 - p_1}{6}$. It gives us a heuristic to find potential answers more quickly.  $p_2$ is in the form of $6k + p_1$.
A: Well, as per your request, I'd add an answer so the comments section is clear (sorry again for the inconvenience).
Notice that you can rearrange the equation as $p = (y + 3x)(y-3x)$.
The least bound for $p$ is $10^9$ as per your question, so I'd always suggest starting an exhaustive search (due to my being a rookie/ newbie at such big problems) from there so that you can get the least possible $p, x$ and $y$ that satisfies your equation. From there, tactically plan how you can manipulate the value of $x$ and $y$ and still obtain values for $p$. For example, you can add a $9a^2 \pm 6ax$ or $a^2 \pm 2ay$ (I'd suggest going for $9a^2 + 18ax$ or $a^2+ 2ay$ as long as you don't override a solution during the exhaustive search) to both sides of the equation and see the changes yourself.
If you're looking for multiple solutions for the same $p$, you can exchange things around.
Edit : OmG's answer is worth a special mention, please go check it above.
Edit: Another great contribution by Servaes; please go check below
A: An experimental approach:
Here I find one solution for this equation which a family of solutions can be based on:
$9x^2=y^2-p=(y-\sqrt p)(y+\sqrt p)$
we can can construct following system of equations:
$\begin{cases}y-\sqrt p=9\\y+\sqrt p=x^2\end{cases}$
Subtracting first equation from second one we get:
$2\sqrt p=x^2-9=(x-3)(x+3)$
Suppose:
$x-3=2\Rightarrow x=5$
$x+3=\sqrt p$
subtracting first fro, second we get:
$\sqrt p-2=6\Rightarrow \sqrt p=8\rightarrow  p=64$
Now we have:
$9x^2=y^2-64=(y-8)(y+8)$
Let $y-8=9\rightarrow y=17$
$x^2=y+8=25$
hence :
$x=5$, $y=17$ and $p=64$
$9\times5^2+64=17^2\space\space\space\space\space (1)$
$10^{12}\div 64=15625\times10^6>10^9$
So we can construct following relation by multiplying both sides of the relation (1) by $15625\times 10^6=(125\times 10^3)^2$:
$9(5\times 125000)^2+15625\times 10^6=(17\times125000)^2$
that is one solution of equation under condition $p=15625000000>10^9$ is:
$x=625\times10^3$
$y=2125\times 10^3$
This experiment shows that only with particular values of p the equation may have integer solutions.
A: If we let $\quad  A^2+B^2=C^2
\qquad A=3x\quad B=\sqrt{p}\quad C=y\qquad$
we can use Euclid's formula for Pythagorean triples to find infinite solutions based on the value of $(x)$.  We start with the formula
$$ \qquad A=m^2-k^2\qquad B=2mk \qquad C=m^2+k^2\qquad$$
and solve the $A$-function for $(k).\quad$
Then we test a defined range of $m$-values to see which yield integers.
\begin{equation}
A=m^2-k^2\implies k=\sqrt{m^2-A}\\
\text{for}\qquad  \lfloor\sqrt{A+1}\rfloor \le m \le \frac{A+1}{2}
\end{equation}
The lower limit ensures $k\in\mathbb{N}$ and the upper limit ensures $m> k.$
There are Pythagorean triples for all $A=2n+1$ but A=$3x$,  so we are limited to odd muliples of $3:$
$$ x\in\big\{1,3,5,7,\cdots\big\}\implies
A\in\big\{3,9,15,21,\cdots\big\}.\quad $$
Here, we will use $x=5\implies A=15$
$$A=15\implies \lfloor\sqrt{15+1}\rfloor=4\le m \le \frac{15+1}{2} =8\\
\land\quad m\in\{4,8\}\implies k \in\{1,7\} $$
$$F(4,1)=(15,8,17)\implies (x,p,y)=(5,8^2,17)\\
 F(8,7)=(15,112,113)\implies (x,p,y)=(5,112^2,113) $$
