Substitution yields the condition on the coefficients $$\left(\frac{\hbar^2}{2m}q^2-\epsilon\right)c_{\mathbf{q}}+\sum_{\mathbf{K'}}U_{\mathbf{K'}}c_{\mathbf{q}-\mathbf{K'}}=0$$
To understand the condition on the $\mathbf{q}$'s you have to make a step back in the derivation of that equation.
Assume
$$\psi(\mathbf{r})=\sum_{\mathbf{q}}c_{\mathbf{q}}e^{i\mathbf{q}\cdot\mathbf{r}},$$
with the $\mathbf{q}$'s arbitrarily chosen. Substituting this ansatz in the Schrödinger equation and rearranging the terms yields
$$\sum_{\mathbf{q}}\left(\frac{\hbar^2}{2m}q^2-\epsilon\right)c_{\mathbf{q}}\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{r}}+\sum_{\mathbf{q},\mathbf{K'}}U_{\mathbf{K'}}c_{\mathbf{q}}\mathrm{e}^{\mathrm{i}(\mathbf{q}+\mathbf{K'})\cdot\mathbf{r}}=0.$$
If the sets $\{\mathbf{q}\}$ and $\{\mathbf{q}+\mathbf{K'}\}$ do not coincide, the exponential functions in the two terms on the left hand side of the above equation are linearly independent, and the only way in which that equation can be satisfied is to have $c_\mathbf{q}=0$, that is, $\psi(\mathbf{r}) = 0$.
If, instead, the sets $\{\mathbf{q}\}$ and $\{\mathbf{q}+\mathbf{K'}\}$ are equal, that is, $\{\mathbf{q}\}$ is invariant for translations of reciprocal lattice vectors, you can relabel $\mathbf{q}+\mathbf{K'}$ as $\mathbf{q}$ and rewrite the second term as
$$\sum_{\mathbf{q},\mathbf{K'}}U_{\mathbf{K'}}c_{\mathbf{q}-\mathbf{K'}}\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{r}},$$
and then collect the exponential $\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{r}}$ between the two terms to obtain
$$\sum_{\mathbf{q}}\left[\left(\frac{\hbar^2}{2m}q^2-\epsilon\right)c_{\mathbf{q}}+\sum_{\mathbf{K'}}U_{\mathbf{K'}}c_{\mathbf{q}-\mathbf{K'}}\right]\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{r}}=0.$$
The above equation holds if and only if the coefficients of the exponentials $\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{r}}$ are zero, and this yields the required condition on the coefficients.
Best Answer
Here's a simple-minded answer:
Let's just compute the momentum of a particle with a Bloch wave function
$$\begin{eqnarray} \left.\langle x \right| \hat{p}\left|\Psi \rangle\right. &=& -i\hbar \left(\frac{d}{dx}\right) u_k(x) e^{i k x} \\ &=& -i \hbar \left( i k u_k(x) e^{ikx} + u_k'(x)e^{ikx}\right) \\ &=& \left( pu_k(x) - i\hbar u_k'(x)\right)e^{ikx} \end{eqnarray}$$
where in the last line we defined $p\equiv \hbar k$. This pretty clearly shows that the Bloch wave function is not an eigenfunction of the momentum operator. So, while you can always break the wave function down into plane waves $e^{ikx}$, and each component is a momentum eigenstate with momentum $p=\hbar k$, the Bloch functions are not themselves momentum eigenstates. Therefore, $k$ in $u_k(x)e^{ikx}$ is not the momentum of the Bloch state. Note, however, that if $u_k(x)=\text{constant}$ so that $u_k'(x)=0$, then we get
$$\left.\langle x \right|\hat{p} \left| \Psi \rangle \right. = pu_k(x)e^{ikx}=p\left.\langle x \ |\Psi \rangle \right. = \left.\langle x \right|\ p \left| \Psi \rangle \right. $$ or in other words
$$\hat{p}\left|\Psi\rangle\right. = p\left|\Psi \rangle\right. .$$
Please make a separate question for the Brillouin zone thing. I would like to answer this, but it belongs in a separate question.