% Here we are going to solve the boxed matrix equation on page one of
% the 1D Weak Period Potential.

% Let K be the smallest reciprocal lattice vector.  All other
% K vectors are integer multiples of K.

K = 1;

% Rescale all energies by epsilon_K = hbar^2 K^2/2m .
% This means that U is actually U/epsilon_K and all energies  
% are in units of epsilon_K.  

% In this simple calculation we take all the U's to be the same
% except for U_0 = 0.

U = 0.5;

% Keep (2*Nmax + 1) reciprocal lattice vectors.
% This can be changed to check for convergence.

Nmax = 5;

Kvalues = (-Nmax:1:Nmax)*K;

% The first Brillouin zone is between -K/2 and K/2.

kvalues = linspace(-K/2,K/2,1000);
energies = zeros(length(Kvalues),length(kvalues));

for counter = 1:length(kvalues);
    k = kvalues(counter);

    H = diag(((k - Kvalues).^2 -U));
    H = H + U*ones(length(Kvalues));

    energy = eig(H);
    energies(:,counter) = energy;
end

plot(kvalues,energies)
axis([min(kvalues) max(kvalues) 0 5])
xlabel('k/K')
ylabel('energy/epsilon_K')
title('U = 0.5, not small perturbation')