% This script is diagonalizing a band matrix
% and showing that the analytic values really
% coincide with the numerical results
% The matrix is 
% [  a  b  0  0  0  ....
%    b  a  b  0  .......
%    0  b  a  b  0  ....
%    0  0  b  a  b  0  .
%    ..... 0  b 
%    .......       a  b
%    0  0  0       b  a ]
%
%%%  eng -> a   V1 -> b
eng=4;V1=-2;
nnn=input('dimension of the matrix: ')  ;
    H0=eng*eye(nnn);
    H1=eye(nnn-1);H1=[zeros(nnn-1,1) H1];H1=[H1;zeros(1,nnn)];H1=H1+H1';
    H1=H0+H1*V1;

[a_vec e_val] = eig(H1);
[vals inds]=sort(diag(e_val));
figure(1)
xx=1:nnn;
plot(xx,vals,'*');hold on
plot(xx,eng+2*V1*cos(pi*xx/(nnn+1)),'-')
ind1=1;
hold off
while (ind1>0) 
   ind1=input('which eigenvector to plot :  (negative stop) ');
   if ind1>0
     vec=a_vec(:,inds(ind1));
     teo=sin(ind1*pi*xx/(nnn+1));
     teo=teo*(vec(1)/teo(1));
     figure(2)
     plot(xx,vec,'o',xx,teo,'-')
   end  
end