clear
n_states=500;
BW=1;
MB=0;
eps0=MB+2.*BW./n_states.*[-n_states./2:1:n_states./2];
H0=diag(eps0);

eps_ad=0.0;
V_ad=2;
H=[eps_ad, V_ad./sqrt(n_states).*ones(1,length(eps0));
    V_ad./sqrt(n_states).*ones(length(eps0),1),H0];
[EV,EE]=eig(H);

figure(2)
plot(diag(EE),EV(1,:).^2,'r.-')

