numerov with  cosine potential

(*Potential*)
V[s_]:=10*Cos[s];
(*fix size of hamiltonian at 500 x 500, and sample s ∈(((-π)/500),2π-(π/500)) uniformly*)
n=500;d=2.π/n;
s=Table[-((d )/(2. ))+d i,{i,n}];
(*calculate KE matrix, impose periodic boundary conditions by coupling ψ[n] and ψ[1] with each other*)
𝕀[n_,d_]:=DiagonalMatrix[1+0Range[n-Abs[d]],d];
B=(1./12)(𝕀[n,-1]+10 𝕀[n,0]+𝕀[n,1]+𝕀[n,n-1]+𝕀[n,-n+1]);
A=(1/( (d^2)))(𝕀[n,-1]-2 𝕀[n,0]+𝕀[n,1]+𝕀[n,n-1]+𝕀[n,-n+1]);
KE=((-1)/2)Inverse[B].A;
(*Hamiltonian*)
H=KE+DiagonalMatrix[V[s]];
(*Energies, wavefunctions*)
{eval,evec}=Eigensystem[H];

numerov[j_]:=Table[{s[[i]],((evec[[j]][[i]])/(d^0.5))},{i,1,n}]

ListPlot[numerov[-15]]

Theoretical Eigenvalues

v=10;

eths=Sort[Table[8^-1*MathieuCharacteristicB[r,4*v]//N,{r,1,50}],Abs[#1]<Abs[#2]&];

ethc=Sort[Table[8^-1*MathieuCharacteristicA[r,4*v]//N,{r,1,50}],Abs[#1]<Abs[#2]&];

(*Theory Eigenvalues*)

eth=Sort[Join[eths,ethc],Abs[#1]<Abs[#2]&];

eth[[;;20]]

(*Numerov Eigenvalues*)

evals=Sort[eval,Abs[#1]<Abs[#2]&];

evals[[;;20]]