You can discritize the Schrödinger form into chain equations, the potential at each
site along the main diagonal V(na) and the hopping parameter or equivalent kinetic energy
along the two main subdiagonals.. Then just diagonlaize the matrix to get all
eigenvalues.(it is symmetric) lambda divided a will be irrational., a the constant spacing of the sites., lambda the period of the potential. You may want to look into
Wannier functions., and the tight binding method,, if you are doubtfull about this proceedure.
Numerically such a spectrum can be seen to be quite intricate, see e.g. my (old) numerical simulations in https://www.researchgate.net/publication/38384527_The_numerical_spectrum_of_a_one-dimensional_Schrodinger_operator_with_two_competing_periodic_potentials
In particular, Bloch waves (which appear for periodic potentials) are completely perturbed: see (possibly) also, https://www.researchgate.net/publication/222569240_Impurity_bands_and_quasi-Bloch_waves_for_a_one-dimensional_model_of_modulated_crystal
Article The numerical spectrum of a one-dimensional Schrödinger oper...
Article Impurity bands and quasi-Bloch waves for a one-dimensional m...