The fundamental solution generates all solutions

First check $(x_n,y_n)$ defined by $x_n+y_n\sqrt d=(x_1+y_1\sqrt d)^n$ is indeed a solution for every $n$: taking conjugates, $x_n-y_n\sqrt d=(x_1-y_1\sqrt d)^n$, so $x_n^2-dy_n^2=(x_n+y_n\sqrt d)(x_n-y_n\sqrt d)=\left[(x_1+y_1\sqrt d)(x_1-y_1\sqrt d)\right]^n=(x_1^2-dy_1^2)^n=1^n=1$.

Now suppose $(x,y)$ is any positive integer solution not of this form; since $x_n\to\infty$ as $n\to\infty$, there is a unique $n$ with $x_n+y_n\sqrt d \le x+y\sqrt d < x_{n+1}+y_{n+1}\sqrt d = (x_n+y_n\sqrt d)(x_1+y_1\sqrt d)$.

Divide through by $(x_n+y_n\sqrt d)$, i.e. multiply by its inverse $(x_n-y_n\sqrt d)$ (valid since $x_n^2-dy_n^2=1$): set $x'+y'\sqrt d = (x+y\sqrt d)(x_n-y_n\sqrt d)$. Then $1\le x'+y'\sqrt d < x_1+y_1\sqrt d$, and $x'^2-dy'^2=(x^2-dy^2)(x_n^2-dy_n^2)=1\cdot 1=1$, so $(x',y')$ is also a solution of Pell's equation.

A short computation using $x'+y'\sqrt d\ge 1$ and $x'^2-dy'^2=1$ shows $x'\ge 1$ and $y'\ge 0$ (a solution with $x'+y'\sqrt d\ge 1$ but $y'<0$ would force $x'>x_1$, contradicting $x'+y'\sqrt d<x_1+y_1\sqrt d$ combined with the norm equation). If $y'>0$, then $(x',y')$ is a positive solution with $x'+y'\sqrt d<x_1+y_1\sqrt d$, contradicting minimality of the fundamental solution. So $y'=0$, forcing $x'=1$, i.e. $x+y\sqrt d=x_n+y_n\sqrt d$, so $(x,y)=(x_n,y_n)$ after all — contradicting our assumption.

Hence every positive solution is exactly some $(x_n,y_n)$, proving the fundamental solution generates the entire solution set.