clear;
hbar2_d_2m=1;
G=[-5 -4 -3 -2 -1 0 1 2 3 4 5]; %in Units of 2.*pi/a_lat
V=[0 0.00 0 0.0 -0.5 0 -0.5 0.0 0 0.00 0];
VMatrix=zeros(length(G));
for g=1:length(G)
    for q=1:length(G)
        g_prime=find(G==(G(g)+G(q)));
        VMatrix(g,g_prime)=VMatrix(g,g_prime)+V(q);
    end
end
%imagesc(VMatrix)
k=-0.5:0.01:0.5 %in Units of 2.*pi/a_lat;

for a=1:length(k);
    H=hbar2_d_2m.*diag((G+k(a)).^2)+VMatrix;
    [EV,EE]=eig(H);
    EVal(a,:)=diag(EE); 
end
k_special=0.0;
H=hbar2_d_2m.*diag((G+k_special).^2)+VMatrix;
[EV,EE]=eig(H);



figure(1)
plot(k,EVal,'*')
ylim([-0.50 3])
figure(2)
x=-8.*pi:0.01:8.*pi;
Pot_RS=0.*x;
WF1=0.*x;
WF2=0.*x;
for a=1:length(G)
Pot_RS=Pot_RS+V(a).*exp(i.*G(a).*x);
WF1=WF1+EV(a,1).*exp(i.*G(a).*x);
WF2=WF2+EV(a,2).*exp(i.*G(a).*x);
end
WF1=WF1./sqrt(sum(abs(WF1).^2));
WF2=WF2./sqrt(sum(abs(WF2).^2));
plot(x./(2.*pi),real(Pot_RS))

title('Potential')

figure(3)
plot(x./(2.*pi),abs(WF1).^2,'r')
hold on
plot(x./(2.*pi),abs(WF2).^2,'g')
hold off
title('Wavefunction')
