Try these monsters:
$x_{n_1}=\frac{1}{2a}\left[\frac{\gamma_0 \sqrt{P^2-1}+Q\beta_0}{2\sqrt{P^2-1}}\cdot \left(P+\sqrt{P^2-1}\right)^n+ \frac{\gamma_0 \sqrt{P^2-1}-Q\beta_0}{2\sqrt{P^2-1}}\cdot \left(P-\sqrt{P^2-1}\right)^n-b\right]$
and
$x_{n_2}=\frac{1}{2a}\left[\frac{\gamma_0 \sqrt{P^2-1}-Q\beta_0}{2\sqrt{P^2-1}}\cdot \left(P+\sqrt{P^2-1}\right)^n+ \frac{\gamma_0 \sqrt{P^2-1}+Q\beta_0}{2\sqrt{P^2-1}}\cdot \left(P-\sqrt{P^2-1}\right)^n-b\right]$
6 notes:
- $(a,b,c)$ from $f(x)=ax^2+bx+c$
- use minimal $(P,Q)$ such that $P^2-aQ^2=1$
- use minimal $(\gamma_0,\beta_0)$ such that $a\gamma_0^2-\beta_0^2=a(b^2-4ac)$
- for non-square $a$, when a is a square it's analogous to the explanation below
- Reasoning for $x_{n_1}$ and $x_{n_2}$ is that some solutions to $ \ a\gamma^2-\beta^2=a(b^2-4ac) \ $ have two separate strands.
- In using this as written, it may be that you need to take every $ \ n=2m \ $ or $ \ n=2m-1 \ $ or multiples $n=km$ to meet the $\gamma-b \equiv 0 \pmod {2a}$ constraint (referring to how $ \ 2a \cdot x + b = \gamma \ $ below)
The derivation of the above is arriving at the pell equation:
$$\begin{align} ax^2 + bx + c &= f(x)=\alpha^2 \\
a(x+b/2a)^2-\frac{b^2-4ac}{4a} &=\alpha^2 \\
a(2ax+b)^2-a(b^2-4ac)&=4a^2\alpha^2=\beta^2 \\
a\gamma^2-\beta^2&=a(b^2-4ac)
\end{align}$$
And solving using standard techniques.Thus
$$x_n=\frac{\gamma_n-b}{2a}$$
Such that $\gamma_n$ is the nth solution in the above pell type equation : $a\gamma^2-\beta^2=a(b^2-4ac)$
An explanation as to why in your first problem $(2,10,19,47)$ are the only answers is what follows:
The problem at the top has $(a,b,c)=(4,84,-15)$. Since $a$ is square, that pell equation morphs into a difference of two squares.
$$ \begin{align}
4\gamma^2-\beta^2&=29184 \\
(2\gamma)^2-\beta^2&=29184=2^9\cdot3\cdot19
\end{align}$$
Using this identity:
$$\left(\frac{d_1+d_2}{2}\right)^2-\left(\frac{d_1-d_2}{2}\right)^2=d_1\cdot d_2$$
Let $d_1, d_2$ be two divisors (of the same parity) of $29184$. You can see here that since there can only be a limited number of divisors of $29184$, even more so limited that have the same parity, you know right here that there's going to be a limited number of solutions.
Let $2\gamma=\frac{d_1+d_2}{2}$, thus $\gamma=\frac{d_1+d_2}{4}$, and finally putting this back into our $x$, and letting $b=84$ and $a=4$, arrive at
$$x=\frac{\left(\frac{d_1+d_2}{4}\right)-84}{8}$$
or
$$x=\frac{d_1+d_2-336}{32}=\frac{d_1+d_2}{32}-\frac{21}{2}$$
Now we observe that both factors must add to something directly proportional to 16:
$$d_1+d_2=16\cdot k$$
for the fraction $21/2$ to become whole when added. This implies that both $d_1$ and $d_2$ need to have a factor of atleast $2^4$. The only two-factor sets that come from $2^9\cdot3\cdot19$ and meet this requirement are:
$$\begin{align}
(d_1,d_2) &\to x(d_1,d_2) &&\to x=\frac{d_1+d_2}{32}-\frac{21}{2} \\
\cdots \cdots \cdots \cdots \cdots \cdots \\
(3\cdot 19 \cdot 2^4, 2^5) &\to x(912,32) &&\to x=19 \\
(3\cdot 19 \cdot 2^5, 2^4) &\to x(1824,16) &&\to x=47 \\
(3\cdot 2^4,19 \cdot 2^5) &\to x(48,608) &&\to x=10 \\
(3\cdot 2^5,19 \cdot 2^4) &\to x(96,304) &&\to x=2 \\
\end{align}$$
Best Answer
As the comments point out, it's not really any faster to use algebraic number theory to solve the problem.
On the other hand, you can consider
$$y^2=-27x^2 + 18px + p^2$$
where $p$ is a prime number (your question is the special case $p=13$). Then solving this equation can be reduced to finding $a$, $b$ such that
$$a^2+ab+b^2=p^2\text{.}$$
Finding solutions to this last equation depends on whether the remainder of $p$ when divided by $3$ is $0$, $1$, or $2$.