Magnitude of band gap on including higher Fourier components

It is commonly shown that the magnitude of the band gap at the edge of the 1st Brillouin Zone is equal to twice the Fourier component of the potential energy for a one dimensional crystal. Also, the way the spectrum is plotted and compared to the free electron dispersion suggests that the split is symmetric. However, all this assumes the underlying potential is a cosine function, and hence a single Fourier component $V_1$.

If I include higher order Fourier components $V_2$ and $V_3$,

1. Will the band 1 and 2 energy gap still be exactly equal to $2V_1$?
2. Will the splitting be symmetric? As in, will the first band get a reduction of energy by $V_1$ and the second band get lifted up by $V_1$ compared to the free electron dispersion?

From my understanding, I expect the band gap at $k=pi/a$ be equal to $2V_1$ and another band gap appear between the third and fourth band of magnitude $2V_3$. Also, at $k=0$ (between 2nd and 3rd band), I would expect a gap magnitude $2V_2$.

What I did:

I expanded the potential and the wavefunction in the plane wave basis and obtained the Hamiltonian matrix for a particular wavevector $k$ in the FBZ. Now, the eigenvalues of the matrix should give me the energies at that particular $k$, which I sort and label them with the band index $n$. Note that the obtained Hamiltonian matrix is a heptadiagonal band matrix since I’ve included the second and third Fourier components. Also, I’ve focussed only on the $5times5$ sub block matrix since I want to investigate the five lowest energies.

I had to resort to numerical methods to find the eigenvalues of the matrix at a particular $k$. I chose the lattice constant of the direct lattice to be $pi$ and set $frac{hbar^2}{2m}=1$. I chose the Fourier components of the potential as $V_0=0, V_1=0.9, V_2=0.4, V_3=0.1$ (in accordance with the fact that the components decrease). Explicitly, my obtained matrix is
$$
begin{pmatrix}
&(k-4)^2 &V_1 &V_2 &V_3 &0
&V_1 &(k-2)^2 &V_1 &V_2 &V_3
&V_2 &V_1 &k^2 &V_1 &V_2
&V_3 &V_2 &V_1 &(k+2)^2 &V_1
&0 &V_3 &V_2 &V_1 &(k+4)^2
end{pmatrix}
$$
Now, for $k=1$ (FBZ boundary), the eigenvalues are obtained using MATLAB. The obtained energies are listed below(1st column) and compared to the free electron energies(2nd column) and what I would have expected(3rd column):

25.0597    25.000    25.0000
 9.2994    9.0000     9.1000
 8.9978    9.0000     8.9000
 1.6683    1.0000     1.9000
 0.0695    1.0000     0.1000

Clearly, the energy eigenvalues are not symmetric and the band gap is not equal to twice the Fourier component. Also, even though the 5th Fourier component of the potential is zero, I get a change in the energy of the 5th band. What is it that I’m doing wrong here?

P.S. I would love a solution to this that doesn’t rely on perturbation theory.