where the periodic function $u_{\vec{k}}=\sum\limits_{\vec{G}'}C_{\vec{G}'}e^{i\vec{G}'\cdot\vec{r}}$ has been expressed as a Fourier series. The reciprocal lattice vectors have been relabeled as running over $\vec{G}'$ instead of $\vec{G}$. It does not matter how they are labeled since we sum over all of the reciprocal lattice vectors, $\sum\limits_{\vec{G}}C_{\vec{G}}e^{i\vec{G}\cdot\vec{r}}=\sum\limits_{\vec{G}'}C_{\vec{G}'}e^{i\vec{G}'\cdot\vec{r}}$. This form for the wavefunction is substituted into the Schrödinger equation,
$$\sum\limits_{\vec{G}'}\dfrac{\hbar^2(\vec{k}+\vec{G}')^2}{2m} C_{\vec{G}'}e^{i(\vec{k}+\vec{G}')\cdot\vec{r}} + \sum_{\vec{G}}\sum\limits_{\vec{G}''} U_{\vec{G}} C_{\vec{G}''}e^{i(\vec{k}+\vec{G}+\vec{G}'')\cdot\vec{r}} = E \sum\limits_{\vec{G}'}C_{\vec{G}'}e^{i(\vec{k}+\vec{G}')\cdot\vec{r}}. $$
In the middle term with the double sum, the sum over $\vec{G}'$ has been relabeled as a sum over $\vec{G}''$. Again it does not matter that the label has changed since the sum is over all of the states but we need a way to keep track of the product terms in the double sum. Next we collect like terms. The exponential factors can be written as $e^{i \vec{k} \cdot \vec{r}}= \cos(\vec{k} \cdot \vec{r})+ i\sin(\vec{k} \cdot \vec{r})$. Only terms that have the same wavelength can be equal to each other so only the terms where $\vec{k} +\vec{G}'= \vec{k} + \vec{G}+\vec{G}''$ can be equal to each other. This results in the condition $\vec{G}'' = \vec{G}' - \vec{G}$ and for each $\vec{G}'$ there is an equation,
$$\dfrac{\hbar^2(\vec{k}+\vec{G}')^2}{2m}C_{\vec{G}'} + \sum_{\vec{G}} U_{\vec{G}} C_{\vec{G}'-\vec{G}} = E C_{\vec{G}'}. $$
This set of equations are called the central equations. The Schrödinger equation, a differential equation for $\psi_{\vec{k}}$, has been replaced with the central equations which are algebraic equations for the coefficients $C_{\vec{G}'}$. The algebraic equations can be put in the form of an eigenvalue problem. These equations written out in the one-dimensional case for $-G_0, 0, G_0$ are,
If the terms indicated by $\cdots$ are negelected, this can be written in matrix form as an eigenvalue problem that can be solved for the energies.
$$\textbf{M} \vec{C} = E \vec{C}.$$
When more equations for $\vec{G}'$ are included, the matrix gets bigger and since the coefficients $U_{\vec{G}}$ typically get smaller as $\vec{G}$ gets larger, it becomes more justified to neglect the terms indicated by $\cdots$. An $N\times N$ matrix will have $N$ solutions for the energy $E$ at each value of $\vec{k}$. The different solutions correspond to different bands.