Which coefficients in the central equation need to be retained for approximate solution at (or near) the zone boundary? To obtain the band gaps near the zone boundaries from the central equation
$$\left(\frac{\hbar^2k^2}{2m}-\varepsilon\right)c_{k}+\sum_G U_G c_{k-G}=0,
$$ where $U_G$ are the Fourier coefficients of the periodic potential $U(x)$.
Why is it adequate to consider only a finite number of $c_k$ coefficients? For example, in Kittel's book, in order to calculate the solution $\varepsilon$ exactly at the zone boundary (i.e. $k=G/2$), only the coefficients $c_{G/2}$ and $c_{-G/2}$ are retained in the central equation, and every other coefficient is neglected. Why? In other words, how do we know that the only important coefficients are $c_{G/2}$ and $c_{-G/2}$?
 A: 
It helps to have the picture above in mind, which shows the free-particle dispersion relation $\epsilon_k = \hbar^2k^2/2m$ plotted in reduced-zone scheme.  Choose a particular $k$ in the first Brillouin zone as shown. Then the central equation shows that only those free-particle states with $k+n2\pi/a$ for some integer $n$ are coupled to each other (indicated by the red dots in the diagram).
Now, the central equation is exact, but it is most useful when we are considering nearly free electrons; that is, the periodic potential $U$ is a perturbation, which means that in general, $U_{G}$ is "small". In that case, according to 2nd-order perturbation theory, the 2nd-order corrections to the energy are given by
$$
E_k^{(2)}=\sum_{n\neq0}
\frac{
\left\lvert
\left\langle\psi_{k+2\pi m/a}^0
\left\lvert \hat{U}\right\rvert
\psi_{k}^0 \right\rangle
\right\rvert^2
}
{\epsilon_k-\epsilon_{k+2\pi m/a}}\,.
$$
(The matrix element is just $U_{2\pi m/a}$.)  This expression tells us that the corrections to the energy decrease with the difference in energy between our "target" state (in this case the free-particle state $\psi_k^0$) and the states that are getting mixed in ($\psi_{k+2\pi m/a}^0$).
We can see from the picture that the state that is closest in energy to the $k$ state (at least on the right side of the first Brillouin zone) is the $k-2\pi/a$ state.  In fact, at the edge of the first Brillouin zone boundary, these states are degenerate, and the other states are very far-removed in energy. This means that these two low-energy states will mix together very strongly while only mixing a little bit of the higher energy states.  Therefore, to a very good approximation, we need only keep those two states when computing the energies of the two lowest-energy states at the Brillouin zone boundary.

More generally, we can perform a variational calculation by constructing a truncated version of the matrix equation implied by the central equation. To do this, we order the basis as shown in the diagram, and construct the matrix whose eigenvalues yield the exact energies of this system. For convenience, let's set $U_0=0$.  We can do this because $U_0$ represents a overall shift of the potential energy, whose absence doesn't affect the physics at all. In addition, we will set $U_G=U_{2\pi m/a}\to U_m$ so that we don't have to carry around a bunch of $2\pi/a$'s. Then,
$$
\begin{bmatrix}
\epsilon_k & U_{-1} & U_{2} & U_{-4} & \cdots \\
U_1 & \epsilon_{k-2\pi/a} & U_{2} & U_{-1} & \cdots \\
U_{-2} & U_{-2} & \epsilon_{k+2\pi/a} & U_{-3} & \cdots \\
U_4 & U_{1} & U_3 & \epsilon_{k-4\pi/a} & \cdots \\
\vdots & \vdots & \vdots & \vdots
\end{bmatrix}
\begin{bmatrix}
c_k\\c_{k-2\pi/a}\\c_{k+2\pi/a}\\c_{k-4\pi/a}\\\vdots
\end{bmatrix}
=
\varepsilon
\begin{bmatrix}
c_k\\c_{k-2\pi/a}\\c_{k+2\pi/a}\\c_{k-4\pi/a}\\\vdots
\end{bmatrix}\,.
$$
By truncating this matrix at some point and diagonalizing the resulting matrix, we are performing a variational calculation with the $c$'s as variational parameters. Choosing a good basis is key for a variational calculation, and the above heuristics generated by second-order perturbation theory tells us which states to use: they are the ones shown in the diagram.
