Is numerical lattice wavefunction smooth? -- graphene tight binding case I tried to follow exactly Sec. II.K [page 112-113, Hamiltonian after Eq. (113)] of the standard Review of Modern Physics paper on graphene, which is a tight-binding model of a graphene stripe under magnetic field.

It's periodic and hence fourier transformed along x, but open along y.
The resulted Landau-level-like energy spectrum looks perfectly fine as in the paper. However, I got confused by the wavefunctions since they look somewhat messy, sawtooth, and not smooth.
I haven't played with tight-binding models much and am not sure if this is correct or not. Probably one don't expect lattice wavefunctions to be smooth at all?
Another question is whether Landau level (LL) degeneracy is in general lifted in lattice models. If so, is the lattice LL a certain superposition of many degenerate LLs, which depends on the lattice model details? 
Here I plot, at a certain $k_x$ around the flat bands, norms of wavefunctions of the lowest 4 armchair bands (from left to right) on one sublattice.

 A: The interaction in almost all tight-binding models has finite range (e. g. only includes nearest-neighbor or next-nearest-neighbor interactions), and therefore the associated matrix-valued trigonometric polynomial is smooth (even analytic). Hence, Bloch functions and energy band functions are locally smooth away from band crossings. 
A: I have one possible explaination. I used the infinite square well Hamiltonian (I didn't pay attention to boundary condititions since it will not matter for the big picture). Then I calculated the eigenvectors/eigenvalues in Mathematica.
n = 50;
d = KroneckerDelta;
H = -Table[d[Abs[i - j] - 1] - 2 d[i - j], {i, 1, 
     n}, {j, 1, n}];
v = Eigenvectors@H;
e = Eigenvalues@H;

The unsorted energy spectrum e looks as follows

Notice two things: firstly the lowest energies are to the right, which tells us that you generally can't trust the order of these energies. You would probably want to sort lowest to highest and sort the eigenvectors along with them. Secondly we would expect a quadratic dependence on $k$ where $k$ is the kth eigenvalue. Here $k$ happens to be the wavenumber like in $\psi=A\sin kx$ as you will see from the plots. So we would expect $E\propto k^2$. You can see that this quadratic dependence holds for low values of $k$ but fails for small wavelengths (large $k$). Again the energies are not sorted so large $k$ is on the left in this plot.
Now let's plot two of these eigenvectors


The first and second picture correspond to $k=2$ and $k=n-2$ respectively. You can see that $k=2$ behaves nicely like we expected. For $k=n-2$ the solution looks like a sawtooth like in your solutions. You can interpret this as the discrete approximation failing for high $k$. When $k$ is high the approximation
$$\frac{d^2\psi}{dx^2}\approx \psi(x+a)-2\psi(x)+\psi(x-a)$$
works less and less well so we can expect the solutions to deviate from the exact answer. I'm not sure if this is what happens in your case but this should be something to keep in mind.
EDIT This probably won't be useful for you anymore given this was asked a year ago but I still answered because this seems to be a problem that more people face.
