Dissatisfied with textbook explanations for why $\vec k$ in Bloch's theorem can be restricted to thefirst Brillouin Zone (BZ)

Consider the 1D Schrödinger's equation with $U(x+a)=U(x)$:

$$-f''(x)+U(x)f(x)-Ef(x)=0,\tag1$$

where $f$ is required (as a boundary condition) to be bounded at infinity: $|f(x)|<\infty$ as $|x|\to\infty$.

Substituting

$$f(x)=u_k(x)\exp(ikx),\tag2$$

we get

$$-u_k''(x)-2iku_k'(x)+(k^2-E+U(x))u_k(x)=0.\tag3$$

By Bloch's theorem, boundary conditions for $u_k$ are

$$ \left\{ \begin{align} u_k(x+a)&=u_k(x),\\ u_k'(x+a)&=u_k'(x). \end{align} \right.\tag4 $$ Now substituting

$$k=k_1+G,\tag5$$

where $G$ is a multiple of $2\pi/a$, we transform $(3)$ into

$$-u_{k_1+G}''(x)-2i(k_1+G)u_{k_1+G}'(x)+((k_1+G)^2-E+U(x))u_{k_1+G}(x)=0.\tag6$$

The boundary conditions for this equation are unchanged, since it's just a replacement of a parameter.

But if now we substitute

$$u_{k_1+G}(x)={u_1}_{k_1}(x)\exp(-iGx),\tag7$$

we'll get an equation

$$-{u_1}_{k_1}''(x)-2ik_1{u_1}_{k_1}'(x)+(k_1^2-E+U(x)){u_1}_{k_1}(x)=0,\tag8$$

which is isomorphic to $(3)$. Moreover, since $\exp(-iGx)$ is periodic with period of $2\pi/a$, boundary conditions are also the same as for $(3)$, i.e. $(4)$. This means that ${u_1}_{k_1}$ and $u_k$ span the same set of solutions.

Now, combining $(2)$, $(5)$ and $(7)$, we get

$$\begin{align} f(x)=u_k(x)\exp(ikx)&=\big[{u_1}_{k_1}(x)\exp(-iGx)\big]\exp(i(k_1+G)x)=\\ &={u_1}_{k_1}(x)\exp(ik_1x), \end{align}\tag9 $$

which is here expressed both in terms of $k$ and in terms of $k_1$, is one and the same solution—for wavenumbers that differ by a whole number of reciprocal lattice constants.