I don't actually know any of the standard proofs of Dirichlet's unit theorem that Qiaochu Yuan refers to, but I think they might use Minkowski's theorem about the existence of nonzero lattice points in sufficiently large centrally symmetric convex subsets of $\mathbf R^n$. I'll give a proof for the specific question asked here, which uses Minkowski's theorem, but only in the very simple case of a parallelogram. I'll write the proof top-down, so that one sees how the theorem is used before I'll state (and prove) it.
First I recall some generalities about the ring $R=\mathbf Z[\sqrt d]$ for $d$ a positive squarefree number, to get started. The additive group of $R$ is free Abelian of rank $2$ with generators $1,\sqrt d$. The norm map $N:R\to\mathbf Z$ given by $a+b\sqrt d\mapsto a^2-db^2$ is multiplicative, and the units of $R$ are precisely the elements with norm $\pm1$. One has the non-trivial unit $-1$, but it is of finite order; the point to prove is therefore that the subgroup of positive units is non-trivial (once a positive unit${}\neq1$ is found, its powers form an infinite set of units). I will reason by contradiction, so assume that $1$ is the unique positive unit of $R$. The first step will be to show that this would imply that for any $n\in\mathbf N$, the number of positive $r\in R$ with $|N(r)|=n$ is finite, in fact at most $n^2$.
Lemma. For any $n\in\mathbf N_{>0}$, the number of principal ideals of $R$ that contain $n$ is at most $n^2$.
Proof. Since these ideals all contain $nR$, they all map to principal ideals of $R/nR$, and the mapping is injective. The number of principal ideals of $R/nR$ cannot exceed the number $n^2$ of its elements. QED
The bound given here is far from sharp, but finiteness is all we need. Under the hypothesis that $1$ is the unique positive unit of $R$, two positive elements of $R$ generate distinct principal ideals, and if $|N(r)|=n$, the ideal generated by $r$ contains $n$, so the lemma justifies our claim that the number of of such$~r$ is finite.
This means that for every $M>0$ there is some $\varepsilon_M>0$ such that for all $a,b\in\mathbf Z$ with $0\neq|a+b\sqrt d|<\varepsilon_M$ one has $|N(a+b\sqrt d)|\geq M$. We shall show that for sufficiently large $M$ this contradicts Minkowski's theorem. Note that $N(a+b\sqrt d)=(a+b\sqrt d)(a-b\sqrt d)$, so we can bound $N(a+b\sqrt d)$ above if in addition to the value $a+b\sqrt d$ we also bound its conjugate $a-b\sqrt d$. Now the linear endomorphism of $\mathbf R^2$ sending $\binom ab\mapsto\binom{a+b\sqrt d}{a-b\sqrt d}$ has determinant $-2\sqrt d$, so the conditions $|a+b\sqrt d|<x$ and $|a-b\sqrt d|<y$ define the interior of a parallelogram of area $4\frac{xy}{2\sqrt d}$, for any $x,y>0$. Minkowski's theorem says such a parallelogram will contain a nonzero lattice point whenever its area is greater than $2^2=4$.
So here is how to obtain a contradiction: take any $M>2\sqrt d$ and put
$$
x_0=\min\{ z\in R \mid z>0 \land |N(z)|\leq M\}\qquad\text{and}\qquad y_0=\frac M{x_0}.
$$
Then $\frac{x_0y_0}{2\sqrt d}>1$, so Minkowski's theorem ensures the existence of $a,b\in\mathbf Z$, not both $0$, with $|a+b\sqrt d|<x_0$ and $|a-b\sqrt d|<y_0$. But then on one hand $|N(a+b\sqrt d)|>M$ by the choice of $x_0$ (and the fact $N(-z)=N(z)$), but on the other hand $|N(a+b\sqrt d)|=|a+b\sqrt d||a-b\sqrt d|<x_0y_0=M$, a contradiction.
Minkowski's theorem. Any centrally symmetric convex subset $S$ of $\mathbf R^d$ of volume greater than $2^d$ contains a nonzero element of $\mathbf Z^d$.
Proof. The natural projection map $f:\mathbf R^d\to(\mathbf R/2\mathbf Z)^d$ is locally volume-preserving, so its restriction to $S$ cannot be injective since the total volume at arrival is $2^d$. If $s,s'\in S$ have $s\neq s'$ and $f(s)=f(s')$ then by central symmetry $-s'\in S$, and by convexity $\frac{s-s'}2\in S$; therefore, since $f(s-s')\in(2\mathbf Z)^d$, one has $\frac{s-s'}2\in S\cap(\mathbf Z^d\setminus\{0\})$. QED
Fermat's equation for cubes is a common introduction to lecture notes on algebraic number theory, because it motivates to study rings of integers in a number field, and partly has been developed even for such Diophantine problems, e.g., Kummer's work concerning generalizing factorization to ideals.
For the equation $x^3+y^3=z^3$ the number field is $\mathbb{Q}(\zeta)$ with a third primitive root of unity $\zeta=e^{2\pi i/3}$. Its ring of integers is given by $\mathbb{Z}[\zeta]$, which is indeed a factorial ring (because it is Euclidean). Its units are given by $\pm 1,\pm \zeta,\pm\zeta^{-1}$. This is crucial to prove Euler's result:
Theorem(Euler $1770$): The equation $x^3+y^3=z^3$ has no non-trivial integer solutions.
The proof uses divisibility properties of the ring $\mathbb{Z}[\zeta]$, starting from the equation
$$
z^3=x^3+y^3=(x+y)(x+\zeta y)(x+\zeta^2y).
$$
The first case is $p=3\nmid xyz$. We may suppose that $x,y,z$ are coprime. We have $z^3\equiv \pm 1\bmod 9$ and $x^3+y^3\equiv -2,0,2 \bmod 9$, so that $x^3+y^3\neq z^3$, a contradiction. Hence we suppose that $3\mid xyz$, i.e., say, $3\mid z$ and $3\nmid xy$. Now we reformulate the equation as
$$
x^3+y^3=(3^mz)^3,
$$
with $x,y,z$ pairwise coprime and $3\nmid xyz$,
where we have solved the case $m=0$. The idea is now to use descent, i.e., to reduce it to the case $m=0$. The above equation becomes
$$
(3^mz)^3=(x+y)(x+\zeta y)(x+\zeta^2y),
$$
where the three factors are not coprime, because $1-\zeta$ is a common factor, because of $3=(1-\zeta)(1-\zeta^2)$, so that $(1-\zeta)\mid 3\mid (x+y)$.
However, since $\mathbb{Z}[\zeta]$ is factorial, all factors are cubes, i.e. ,
$x+y=3^{3m-1}c^3$ with some $c\in \mathbb{Z}$, and so on. This finishes the proof, after some computations in this ring.
Unfortunately, this idea does not work for $x^p+y^p=z^p$ for primes $p$, except for $p\le 19$, because otherwise the ring of integers $\mathbb{Z}[\zeta_p]$ is no longer factorial.
Best Answer
Let's take a solution $(x,y)$ of $x^2-5y^2=\pm4$. Assume $x>0$ and $y>0$. Clearly also $x$ and $y$ have the same parity. Define $$x'=\frac{5y-x}2,\qquad y'=\frac{x-y}2.$$ Then $x'$ and $y'$ are integers, and $$x'^2-5y'^2=\frac{(5y-x)^2-5(x-y)^2}4=\frac{20y^2-4x^2}4=\pm4.$$ Therefore $(x',y')$ is also a solution. I claim that $y'\ge0$ and $x'>0$.
If $y'<0$ then $x<y$ and $x^2-5y^2<-4x^2<-4$, which is false. So $y'\ge0$. If $x'\le0$ then $x\ge5y$ and $x^2-5y^2\ge20y^2>4$ which is false. So $y'\ge0$.
I claim that as long as $y\ge2$, then $y'<y$. Otherwise, $x\ge3y$ and $\pm 4=x^2-5y^2\ge4y^2$. This is only possible if $y=1$.
So, iterating the operation $(x,y)\mapsto(x',y')$ eventually reduces $(x,y)$ to a solution $(X,Y)$ with $X>0$ and $Y\in\{0,1\}$. Therefore to $(X,Y)=(2,0)$, $(1,1)$ or $(3,1)$. All of these reduce down to $(2,0)$.
Therefore we can start with $(x_0,y_0)=(2,0)$ and reversing the operation generate all positive solutions. The iterative process is $$(x_{n+1},y_{n+1})=\left(\frac{x+5y}2,\frac{x+y}2\right).$$