William LeVeque's Topics in Number Theory, Volume I treats Pell's equation with nonsquare $d$ in Chapter 8, without the assumption that $d$ is squarefree.
In Section 8.2 (page 139 in the first half of the Dover edition), he considers $x^2-dy^2 = 1$ as follows:
If $d$ is negative, then when $d=-1$ the only solutions are $(\pm 1,0)$ and $(0,\pm 1)$. If $d\lt -1$, then the only solutions are $(\pm 1,0)$.
If $d$ is a square, d=(d')^2, then we can rewrite the equation as x^2 - (d'y)^2 = 1; the only squares that differ by $1$ are $0$ and $1$, so again the only solutions are $(\pm 1,0)$.
For positive nonsquare $d$, we start with a lemma on approximation and then proceed.
Lemma. If $\xi$ is a real number and $t$ is a positive integer, then there are integers $x$ and $y$ such that $\left|\xi - \frac{x}{y}\right| \leq \frac{1}{y(t+1)},\qquad 1\leq y\leq t.$
Proof. The $t+1$ numbers $0\cdot \xi - \lfloor 0\cdot \xi\rfloor,\quad 1\cdot\xi - \lfloor 1\cdot \xi\rfloor, \quad\cdots\quad t\xi - \lfloor t\xi\rfloor$ are all in the interval $0\leq u \lt 1$. In increasing order of magnitude, call them $a_0$, $a_1,\ldots,a_t$. Mark the numbers on a circle of perimeter $1$. Then the $t+1$ differences $a_1-a_0,\quad a_2-a_1,\quad\ldots,\quad a_t-a_{t-1},\quad a_0-a_t+1,$ are the lengths of the arcs between successive $a$s, so they are nonnegative, and $(a_1-0) + (a_2-a_1) + \cdots + (1-a_t) = 1.$ So at least one of these $t+1$ differences is no more than $\frac{1}{t+1}$. But the each difference is of the form $g_1\xi - g_2\xi - m$, with $m$ an integer, so take $y=|g_1-g_2|$, $x=\pm m$. $\Box$
Theorem. For any irrational $\xi$, the inequality $|x - \xi y| \lt \frac{1}{y}$ has infinitely many integer solutions.
Proof. For each positive integer $t$, $0\lt |x-\xi y|\lt \frac{1}{t}$, $1\leq y\leq t$ has solutions. Taking $t=1$ gives a solution $(x_1,y_1)$ to the original equation. For sufficiently large $t_1$ we have $|x_1 - \xi y_1|\gt \frac{1}{t_1}$, so taking $t=t_1$ gives a new solution $(x_2,y_2)$ to the original equation. Lather, rinse, repeat. $\Box$
The irrationality of $\xi$ is used to ensure that you have strict inequalities when you apply the lemma, so that you can set up the recursion and get infinitely many solutions.
Theorem. If $d$ is positive and not a square, then there are infinitely many integer solutions to the equation $x^2 - dy^2 = k$ for positive integers $x,y$ for some $k$ with $|k|\lt 1+2\sqrt{d}$.
Proof. Pick a solution to $|x - \sqrt{d}y|\lt \frac{1}{y}$. Then \begin{align*} |x+y\sqrt{d}| &= |x-y\sqrt{d} + 2y\sqrt{d}|\\ &\lt \frac{1}{y}+2y\sqrt{d}\\ &\leq (1+2\sqrt{d})y \end{align*} and so $|x^2-dy^2| \lt \frac{1}{y}(1+2\sqrt{d})y = 1+2\sqrt{d}$.
Since there are infinitely many pairs $(x,y)$ you can use, but only finitely many integers that are positive and smaller than $1+2\sqrt{d}$, infinitely many of the values of $x^2-dy^2$ must coincide, which is the theorem. $\Box$
Notice that the only requirement here is that $\sqrt{d}$ be real and irrational, i.e., that $d$ is positive and not a perfect square.
Theorem. If $d\gt 0$ is not a square, then the equation $x^2 - dy^2 = 1$ has at least one solution with $y\neq 0$.
Proof. Take the infinitely many solutions to $x^2-dy^2 = k$ (for some $k$ with $|k|\lt 1+2\sqrt{d}$) and divide them into $k^2$ equivalence classes, where $(x_1,y_1)\sim (x_2,y_2)$ if and only if $x_1\equiv x_2\pmod{k}$ and $y_1\equiv y_2\pmod{k}$. Some class contains more than one solution, say $(x_1,y_1)$ and $(x_2,y_2)$ with $x_1x_2\gt 0$. Let $x = \frac{x_1x_2 - dy_1y_2}{k},\qquad y=\frac{x_1y_2 - x_2y_1}{k};$ then $x$ and $y$ are integers with $y\neq 0$ and $x^2-dy^2 = 1$.
Indeed, $x_1y_2\equiv x_2y_1\pmod{k}$, so $y$ is an integer. Also, $x_1x_2 - dy_1y_2 \equiv x_1^2 - dy_1^2 = k \equiv 0 \pmod{k}$ so $x$ is an integer. Also \begin{align*} x^2 - dy^2 &= \frac{1}{k^2}\left( (x_1x_2 - dy_1y_2)^2 - d(x_1y_2-x_2y_1)^2\right)\\ &= \frac{1}{k^2}\left(x_1^2x_2^2 - dx_1^2y_2^2 + d^2y_1^2y_2^2 - dx_2^2y_1^2\right)\\ &= \frac{1}{k^2}\left(x_1^2-dy_1^2\right)\left(x_2^2 - dy_2^2\right) = 1. \end{align*} Finally, if $y=0$, then $x_1y_2=x_2y_1$, so $x_1=ax_2$ and $y_1=ay_2$ for some $a$; but plugging into $x^2- dy^2 = k$ we get $a=1$, contradicting that $(x_1,y_1)$ and $(x_2,y_2)$ are distinct solutions. $\Box$
Theorem. If $(x_1,y_1)$ and $(x_2,y_2)$ are solutions to $x^2-dy^2 = 1$, then so are the integers $x$ and $y$ defined by $(x_1 + y_1\sqrt{d})(x_2+y_2\sqrt{d}) = x+y\sqrt{d}.$