  We study the Landau level spectrum of bulk graphene monolayers beyond the
Dirac Hamiltonian with linear dispersion. We consider an effective Wannier-like
tight-binding model obtained from ab initio calculations, that includes
long-range electronic hopping integral terms. We employ the Haydock-Heine-Kelly
recursive method to numerically compute the Landau level spectrum of bulk
graphene in the quantum Hall regime and demonstrate that this method is both
accurate and computationally much faster than the standard numerical approaches
used for this kind of study. The Landau level energies are also obtained
analytically for an effective Hamiltonian that accounts for up to third nearest
neighbor hopping processes. We find an excellent agreement between both
approaches. We also study the effect of disorder on the electronic spectrum.
Our analysis helps to elucidate the discrepancy between theory and experiment
for the high-energy Landau levels energies.
