The main book you want is Binary Quadratic Forms: Classical Theory and Modern Computations by Duncan A. Buell. He does not actually quote Lagrange's theorem on bounds, for that Introduction to the Theory of Numbers by Leonard E. Dickson. On page 111 of Dickson, we have Theorem 85 due to Lagrange, any primitively represented value $n$ occurs as the first coefficient of a reduced form in the cycle as long as $|n| < \frac{1}{2} \sqrt \Delta.$
From your edit after Gerry's last suggestion, it seems you may expect to figure this out without looking at any books. Although I give a complete algorithm below, you still will need to do some reading.
Please see my Numbers representable as $x^2 + 2y^2$ , I do not wish to retype everything.
Note: to actually find $\sqrt 8 \pmod p,$ see TONELLI or that other one
CIPOLLA. Or Berlekamp's factoring modulo primes, applied to $w^2 - 8.$
Other note: for an odd prime $q$ such that $(8|q)= -1,$ then whenever $x^2 - 8 y^2 \equiv 0 \pmod q,$ it follows that $x,y \equiv 0 \pmod q,$ so $x^2 - 8 y^2 \equiv 0 \pmod {q^2},$ and $(x/q)^2 - 8 (y/q)^2 $ is a smaller integer (in absolute value anyway).
Every odd prime $p$ such that $(8|p )= 1$ can be represented by either $x^2 - 8 y^2$ or $8 x^2 - y^2.$ These are in distinct genera.
If $p \equiv 1 \pmod 8, $ we can solve $w^2 \equiv 8 \pmod p$ and $p t - w^2 = -8.$ The indefinite quadratic form $\langle p, 2 w, t \rangle$ can be reduced to $x^2 - 8 y^2$ by a matrix in $SL_2 \mathbb Z.$ The inverse of that matrix tells us how to write one solution to $\delta^2 - 8 \gamma^2 = p.$ Then
$$ (\delta^2 + 8 \gamma^2)^2 - 8 (2 \delta \gamma)^2 = p^2. $$
If $p \equiv 7 \pmod 8, $ we can solve $w^2 \equiv 8 \pmod p$ and $p t - w^2 = -8.$ The indefinite quadratic form $\langle p, 2 w, t \rangle$ can be reduced to $8x^2 - y^2$ by a matrix in $SL_2 \mathbb Z.$ The inverese of that matirx tells us how to write one solution to $8 \delta^2 - \gamma^2 = p,$ or $ \gamma^2 - 8 \delta^2 = -p.$ Then
$$ (\gamma^2 + 8 \delta^2)^2 - 8 (2 \delta \gamma)^2 = (-p)^2 = p^2. $$
Given any solution $x^2 - 8 y^2 = q, $ we get a new solution by negating either $x$ or $y.$ Other then that, infinitely many new solutions are created by multiplying by the automorph as in
$$
\left( \begin{array}{cc}
3 & 8 \\
1 & 3
\end{array}
\right)
\left( \begin{array}{c}
x \\
y
\end{array}
\right).
$$
Written as a horizontal formula
$$ (3 x + 8 y)^2 - 8 (x + 3 y)^2 = x^2 - 8 y^2. $$
You can then multiply by the automorph again and again. Note that, because there is no $xy$ term, we can be satisfied with negating entries as we like and do not need to introduce the inverse of the automorph, but it does no harm if you also use it.
I will need to think a bit about completeness for this process. But experiment with it for a while, see what you can do.
POSTSCRIPT: Note that $\langle 1,0,-8 \rangle$ is not actually reduced in the sense of Gauss et al. However, it is one simple step away from $\langle 1,4,-4 \rangle,$ which is reduced.
=-=-=-=-=-=-=
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$ ./indefCycle
Input three coefficients a b c for indef f(x,y)= a x^2 + b x y + c y^2
1 0 -8
0 form 1 0 -8 delta 0
1 form -8 0 1 delta 2
2 form 1 4 -4
-1 -2
0 -1
To Return
-1 2
0 -1
0 form 1 4 -4 delta -1
1 form -4 4 1 delta 4
2 form 1 4 -4
minimum was 1rep 1 0 disc 32 dSqrt 5.6568542495 M_Ratio 32
Automorph, written on right of Gram matrix:
-1 -4
-1 -5
Trace: -6 gcd(a21, a22 - a11, a12) : 1
=========================================
jagy@phobeusjunior:~/old drive/home/jagy/Cplusplus$
=-=-=-=-=-=-=
Let me paraphrase your question as follows:
Determine all $z\in\Bbb{Z}$ for which there exist $d,x,y\in\Bbb{Z}$ with $|d|,|y|>1$ and $d$ squarefree such that
$$x^2+dy^2=z.\tag{1}$$
First note that for $z=0$ there are no integral solutions with $d$ squarefree.
If $z\neq0$ then for every integer $x$ we have the trivial solution
$$(d,x,y)=(z-x^2,x^2,1),$$
which of course fails to meet the condition that $|y|>1$. But for sufficiently large values of $x$ we get
$$d=z-x^2<-1,$$
and so $\Bbb{Z}[\sqrt{-d}]$ is a real quadratic ring. By Dirichlet's unit theorem its unit group has rank $1$, so if $u+v\sqrt{-d}\in\Bbb{Z}[\sqrt{-d}]$ is a fundamental unit and $n\in\Bbb{Z}$ is any integer we have
$$N\left((x+y\sqrt{-d})(u+v\sqrt{-d})^n\right)=N(x+y\sqrt{-d})N(u+v\sqrt{-d})^n=z,$$
yielding infinitely many integral solutions to $(1)$: If $a_n,b_n\in\Bbb{Z}$ are such that
$$a_n+b_n\sqrt{-d}=(x+y\sqrt{-d})(u+v\sqrt{-d})^n,$$
then the above shows that indeed
$$a_n^2+db_n^2=z.$$
Moreover, this yields infinitely many integral solutions $(d,x,y)=(d,a_n,b_n)$ with $|y|>1$, because if $b_m=b_n$ then it quickly follows that $m=n$.
All that remains to be shown is that we can choose $x$ sufficiently large such that $d=z-x^2$ is squarefree.
Best Answer
This is called the generalized Pell equation. As in the classical case there is an algorithm, based on simple continued fractions, due to Lagrange, which solves $$ x^2-dy^2=n $$ for any given squarefree $d$ and given $n\in \Bbb Z\setminus \{0\}$.
Reference: Section $6$ of Keith Conrad's notes.
I would not call this algorithm "trivial" but certainly it is well-known and easy to perform. For small $d$, like $d=2$ it might be a bit quicker, but still is non-trivial.