Here is the sort of thing that is not immediately visible in the quadratic field viewpoint. Suppose we have target $n,$ say $\gcd(n,4d) = 1$ and suppose we are able to solve
$$ u^2 \equiv 4 d \pmod {4n}, $$
$$ u^2 = 4d + 4 n t. $$ Then
$$ u^2 - 4 n t = 4 d. $$ Thus we have constructed a quadratic form
$$ \langle n,u,t \rangle $$
of discriminant $4d.$ I meant to guarantee that the form was primitive. By the reduction scheme of Gauss and Lagrange, this form reduces over $SL_2 \mathbb Z$ to some "reduced" form $ \langle a,b,c \rangle $ of the same discriminant, where reduced is equivalent to $ac <0, b > |a+c|.$ This need not be the principal form, and if $n$ is prime there is just one possible form and its "opposite" that represent $n.$
Here is a good one: we can represent $x^2 - 229 y^2 = p$ (for $p \neq 2,3, 229$) if and only if $z^3 - 4 z - 1 $ factors into three distinct factors $\pmod p.$ This is from an Henri Cohen book, appendix. The form class group, reduced, are
916 factored 2^2 * 229
1. 1 30 -4 cycle length 10
2. 3 28 -11 cycle length 18
3. 11 28 -3 cycle length 18
form class number is 3
=====================================================================
Small positive primes x^2 - 229 y^2 :
37 53 173 193 229 241 347 359 383 439
443 449 461 503 509 541 593 607 617 619
643 691 907 967 977
parisize = 4000000, primelimit = 500509
? factormod( z^3 - 4 * z - 1, 37)
%1 =
[Mod(1, 37)*z + Mod(8, 37) 1]
[Mod(1, 37)*z + Mod(13, 37) 1]
[Mod(1, 37)*z + Mod(16, 37) 1]
? factormod( z^3 - 4 * z - 1, 53)
%2 =
[Mod(1, 53)*z + Mod(5, 53) 1]
[Mod(1, 53)*z + Mod(19, 53) 1]
[Mod(1, 53)*z + Mod(29, 53) 1]
? factormod( z^3 - 4 * z - 1, 173)
%3 =
[Mod(1, 173)*z + Mod(26, 173) 1]
[Mod(1, 173)*z + Mod(156, 173) 1]
[Mod(1, 173)*z + Mod(164, 173) 1]
?
============================================
Monday, 19 Feb.: now i remember why I stuck with positive binary forms for my inhomogeneous examples; I just ran $$3 x^2 + 13 x y - 5 y^2 + z^3 - 4 z \neq \pm 1$$ which is very nice. However, in order to have the discriminant of $z^3 - 4 z + r$ to be a square multiple of that for $r=1,$ we get $27 r^2 + 229 w^2 = 256,$ so we get no more. Compare https://mathoverflow.net/questions/12486/integers-not-represented-by-2-x2-x-y-3-y2-z3-z and Kevin Buzzard's answer. The collection of related material, as pdfs, is now at ZAKUSKI INHOM
Probably enough. I like Buell, Binary Quadratic Forms. Similar in Dickson, Introduction to the Theory of Numbers (1929). Comparison of definitions of "reduced" is in Franz Lemmermeyer's 2010 book, page 37 (pdf 43) Theorem 1.36.
For actually finding all representations of some $n$ of intermediate size by an indefinite form, i like Conway's topograph method. It displays both the action of the (oriented) automorphism group of the form and the finite set of representations that are distinct under the group action. Here is one with pretty good illustrations, I found that drawing by hand on paper works best by drawing the "river" on one page (if it will fit) but then draw "trees" leaving the riverside to show the targets http://math.stackexchange.com/questions/739752/how-to-solve-binary-form-ax2bxycy2-m-for-integer-and-rational-x-y/739765#739765
http://www.maths.ed.ac.uk/~aar/papers/conwaysens.pdf and
https://www.math.cornell.edu/~hatcher/TN/TNbook.pdf and
http://www.springer.com/us/book/9780387955872
http://bookstore.ams.org/mbk-105/
Let $x_n+y_n\sqrt D=(x_1+y_1\sqrt D)^n$ the $n^{\text{th}}$ power of the primitive unit. Since there are only $41^2=1681$ possibilities for $(x_n,y_n)$ $\pmod{41}$ a duplicate must be encountered at some point: $x_n\equiv x_m\pmod{41}$ and $y_n\equiv y_m\pmod{41}$ for some $n>m\ge1$. Then $x_{n-m}=x_nx_m-Dy_ny_m\equiv x_n^2-Dy_n^2\equiv1\pmod{41}$ and $y_{n-m}=-x_ny_m+y_nx_m\equiv-x_ny_n+y_nx_n\equiv0\pmod{41}$.
EDIT: As an example let $D=3$ and the first solution to Pell's equation is $x_1+y_1\sqrt D=2+1\sqrt3$. Now let's make a table of values $\pmod{41}$:
$$\begin{array}{r|r|r}n&x_n&y_n\\\hline
1&2&1\\
2&7&4\\
3&26&15\\
4&15&15\\
5&34&4\\
6&39&1\\
7&40&0\\
8&39&40\\
9&34&37\\
10&15&26\\
11&26&26\\
12&7&37\\
13&2&40\\
14&1&0\\
15&2&1\end{array}$$
For example $(2+1\sqrt3)^2=7+4\sqrt3$, $(2+1\sqrt3)^3=26+15\sqrt3$, and $(2+1\sqrt3)^4=97+56\sqrt3$ so $x_4=97\equiv15\pmod{41}$ and $y_4=56\equiv15\pmod{41}$, thus explaining the row $n=4$, $x_n\equiv15$, $y_n\equiv15$. The first duplicate was $x_{15}\equiv x_1\equiv2\pmod{41}$ and $y_{15}\equiv y_1\equiv1\pmod{41}$, so that tells us that $x_{15-1}=x_{14}\equiv1\pmod{41}$ and $y_{15-1}=y_{14}\equiv0\pmod{41}$. Perhaps a bit anticlimactic since we already found $2$ solutions on our way to generating the first duplicate. Indeed $x_{14}^2-3y_{14}^2=50843527^2-3\cdot29354524^2=1$ and $y_{14}=29354524=41\cdot715964$.
EDIT: Oh yeah, the last $2$ lines: since $(x_n+y_n\sqrt D)(x_n-y_n\sqrt D)=(x_1+y_1\sqrt D)^n(x_1-y_1\sqrt D)^n=(x_1^2-Dy_1^2)^n=(1)^n=1$ we see that $(x_n+y_n\sqrt D)^{-1}=(x_n-y_n\sqrt D)$ so $(x_{n-m}+y_{n-m}\sqrt D)=(x_n+y_n\sqrt D)(x_m-y_m\sqrt D)=(x_nx_m-Dy_ny_m)+(-x_ny_m+y_nx_m)\sqrt D$
EDIT My program that finds the fundamental solution to $x^2-Dy^2=1$ and the first power $n-m$ for which $x_{n-m}\equiv1\pmod{41}$ and $y_{n-m}\equiv0\pmod{41}$
program pell
use ISO_FORTRAN_ENV
implicit none
integer(INT64) D
integer(INT64) sqD, r, s, a, p0, p1, p, q0, q1, q, n
integer(INT64) m
write(*,'(a)') ' D x_1 y_1 n-m'
do D = 1, 100
sqD = int(sqrt(D+0.5D0),INT64)
if(sqD**2==D) cycle
r = 0
s = 1
p0 = 0
p1 = 1
q0 = 1
q1 = 0
do n = 1, 200
a = (sqD+r)/s
p = a*p1+p0
p0 = p1
p1 = p
q = a*q1+q0
q0 = q1
q1 = q
r = a*s-r
s = (D-r**2)/s
if(mod(n,2) == 0 .AND. s == 1) then
write(*,'(i4,1x,i17,1x,i18)',advance='no') D,p,q
p0 = mod(p,41)
q0 = mod(q,41)
p1 = 1
q1 = 0
do m = 1, 1000000
p = p1*p0+D*q1*q0
q = p1*q0+q1*p0
p1 = mod(p,41)
q1 = mod(q,41)
if(p1 == 1 .AND. q1 ==0) then
write(*,'(1x,i4)') m
exit
end if
end do
exit
end if
end do
end do
end program pell
And its output:
D x_1 y_1 n-m
2 3 2 5
3 2 1 14
5 9 4 20
6 5 2 42
7 8 3 21
8 3 1 5
10 19 6 20
11 10 3 42
12 7 2 7
13 649 180 14
14 15 4 7
15 4 1 21
17 33 8 42
18 17 4 5
19 170 39 42
20 9 2 20
21 55 12 40
22 197 42 42
23 24 5 10
24 5 1 42
26 51 10 42
27 26 5 14
28 127 24 21
29 9801 1820 14
30 11 2 42
31 1520 273 5
32 17 3 5
33 23 4 40
34 35 6 21
35 6 1 42
37 73 12 20
38 37 6 42
39 25 4 40
40 19 3 20
41 2049 320 82
42 13 2 40
43 3482 531 10
44 199 30 21
45 161 24 10
46 24335 3588 20
47 48 7 7
48 7 1 7
50 99 14 5
51 50 7 20
52 649 90 14
53 66249 9100 14
54 485 66 14
55 89 12 7
56 15 2 7
57 151 20 40
58 19603 2574 42
59 530 69 10
60 31 4 21
61 1766319049 226153980 5
62 63 8 20
63 8 1 21
65 129 16 42
66 65 8 10
67 48842 5967 42
68 33 4 42
69 7775 936 14
70 251 30 42
71 3480 413 21
72 17 2 5
73 2281249 267000 20
74 3699 430 20
75 26 3 14
76 57799 6630 21
77 351 40 40
78 53 6 8
79 80 9 7
80 9 1 20
82 163 18 82
83 82 9 4
84 55 6 40
85 285769 30996 2
86 10405 1122 20
87 28 3 40
88 197 21 42
89 500001 53000 42
90 19 2 20
91 1574 165 40
92 1151 120 5
93 12151 1260 7
94 2143295 221064 3
95 39 4 7
96 49 5 21
97 62809633 6377352 42
98 99 10 5
99 10 1 42
Best Answer
If $d-1$ is divisible by $4$, an integral element not in $\mathbb{Z}[\sqrt{d}]$ is some $\omega=\frac{a+b\sqrt{d}}{2}$ with $a,b$ odd. Its norm is the integer $N=\frac{a^2-db^2}{4}$.
Then $\omega^2=-N+a\omega$, and $\omega^3=-N\omega+a\omega^2=-aN+(a^2-N)\omega$.
If $N$ is odd, then indeed $\omega^3 \in \mathbb{Z}[\sqrt{d}]$.
Note that for instance $\frac{1+\sqrt{5}}{2}$ is a unit in $\mathcal{O}_{\mathbb{Q}(\sqrt{5})}$ of norm $-1$, and its square isn’t in $\mathbb{Z}[\sqrt{5}]$.