Ashcroft and Mermin are occasionally not too careful about notation, and so the resolution is, as Marek mentioned in a comment, related to the fact that
$$ u_k = \sum_G c_{k-G}e^{iGr} $$
admits several solutions -- every set of coefficients $\{c\}$ corresponds to some function $u_k$, but we're only interested in the $n$ solutions that solve the Schroedinger equation.
Every delocalized band arises because of the hybridization of quantum states between unit cells. If there are $n$ states that the electron can occupy, their will be $n$ (possibly degenerate, depending on the symmetry of the $G$-space point) bands. I say this to clarify that there is no sense in which the band index corresponds to the set reciprocal lattice vectors.
However, to answer your last question, you're absolutely right that reduction to the first Brillouin zone does give you multiple values of $E$ for a given $k \in IBZ$. If I remember correctly A&M introduce band theory by starting with the nearly free electron model and should show explicitly that this is the case. In the tight binding picture (sometimes called Linear Combination of Atomic Orbitals) the energies for $k$ outside the IBZ fall back onto their values in the IBZ.
Continued:
Specifically, I wanted to point out that $E_n(k)=E(k+G_n)$ simply doesn't make any sense. Where did you read that? However the point is that the eigenvalues that make up a particular band are periodic in $G$, so it is true that $E_n(k)=E_n(k+G)$. Now, what this means is that a description with any more than the first Brillouin zone is redundant. The fact that you have multiple bands does not come from the reduction, though it can look that way at first, as A&M present it. I suggest skipping ahead to the chapter on tight binding to understand how multiple bands arise.
Continued again:
Okay, I didn't address this as clearly as I'd hoped. The reason a new quantum number appears at all is that we know from experience that $E(k)$ is not single-valued. The reason for this is that each atom has a set of discrete orbitals, and mathematically we can describe the motion of an electron in a lattice as tunnelings between these orbitals... However, these orbitals have very different symmetries, and we can not, for example, hop from an s-orbital on site-1 to a p-orbital on site-2. Thus each flavor of orbital symmetry will basically correspond to a new band (with a lot of technical complications that will only muddy the point).
Lastly, you're right about $E_n$ being bounded everywhere, and your reasoning is right as well. What was the question?
Right, so it turns out this is quite simple. See the figure below.
(a) shows the 1st BZ resulting from the hexagonal and rectangular unit cells of the honeycomb lattice. Note that the hexagon has twice the area of the rectangle, so two bands from the rectangular BZ will be combined to form one band in the hexagonal BZ.
(b1,2) are a pair of bands calculated for the rectangular unit cell. (Calculate them any way you see fit. For this example, I used a finite difference method and zero potential.) It's easy to identify them as a pair because they have the same values around the edges of the rectangular BZ.
(c) cut the second band into four parts as indicated. In this case, it's pretty clear where to cut.
(d) stitch everything together as indicated. The first band from the rectangular BZ goes in the middle and the cut-up parts surround it.
Done!
Best Answer
This is how you do it: draw your crystal, pick one atom in the crystal and write down the coulomb energy for that atom by considering all neighbouring atoms (distances, charges and the number). See attached manuscript.
Edit: I have realised that I have made a small mistake in the manuscript (the last $4/\sqrt{5}$ in the parentheses needs to be $8/\sqrt{5}$) but to show the principle I kept it.
I have attached a better drawing below to see the exact nearest neighbour distances and count the number of neighbours. The term in parentheses (the Madelung Constant) is an infinite series expansion. You should find out which series is that. For example if you calculate Madelung constant is one dimension the expansion is for $$ln(1+x)= x-x^2/2+x^3/3-x^4/4+...$$