a = 1;
U = 20.0;

[kx,ky] = meshgrid(pi*(-1:0.05:1)/a);

K1x = 2*pi/a; K1y = 2*pi/a;
K2x = 2*pi/a; K2y =-2*pi/a;
K3x =-2*pi/a; K3y =-2*pi/a;
K4x =-2*pi/a; K4y = 2*pi/a;

Npts = size(kx,1); % Assumes square in k-space
band1 = 0*kx;
band2 = 0*kx;
band3 = 0*kx;
band4 = 0*kx;
band5 = 0*kx;

for mx = 1:Npts,
    for my = 1:Npts,
        qx = kx(mx,my);
        qy = ky(mx,my);
        ep0q = qx^2 + qy^2;
        ep0qK1 = (qx - K1x)^2 + (qy - K1y)^2;
        ep0qK2 = (qx - K2x)^2 + (qy - K2y)^2;
        ep0qK3 = (qx - K3x)^2 + (qy - K3y)^2;
        ep0qK4 = (qx - K4x)^2 + (qy - K4y)^2;
        H = [ ep0q U/4 U/4 U/4 U/4;
              U/4 ep0qK1 0 0 0;
              U/4 0 ep0qK2 0 0;
              U/4 0 0 ep0qK3 0;
              U/4 0 0 0 ep0qK4 ];
        E = sort(eig(H)); % sorting places lowest at 1 position
        band1(mx,my) = E(1);
        band2(mx,my) = E(2);
        band3(mx,my) = E(3);
        band4(mx,my) = E(4);
        band5(mx,my) = E(5);
    end
end

figure(1)
hold off
surf(kx,ky,band1)
hold on
surf(kx,ky,band2)
surf(kx,ky,band3)
hold off

figure(2)
contour(kx,ky,band1)
xlabel('kx')
ylabel('ky')
title('band 1')

figure(3)
contour(kx,ky,band2)
xlabel('kx')
ylabel('ky')
title('band 2')
