Nmax = 2;
U = 0.1;
[Kx,Ky] = meshgrid(-Nmax:1:Nmax);
kvalues = linspace(-1/2,1/2,100);
[kxvalues kyvalues] = meshgrid(kvalues);

band1 = zeros(length(kvalues));
band2 = zeros(length(kvalues));
band3 = zeros(length(kvalues));

for cx = 1:length(kxvalues);
    for cy = 1:length(kyvalues);

    kx = kxvalues(cx,cy);
    ky = kyvalues(cx,cy);

    epsilonk0 = (kx - Kx).^2 + (ky - Ky).^2;
    epsilonk0 = epsilonk0(:);

    H = diag(epsilonk0 -U);
    H = H + U*(0*H + 1);

    energy = eig(H);
    energy = sort(energy);
    band1(cx,cy) = energy(1);
    band2(cx,cy) = energy(2);
    band3(cx,cy) = energy(3);
    end
end

figure(1)
contour(kxvalues,kyvalues,band1)
xlabel('k_x')
ylabel('k_y')
title('band 1')
axis('square')

figure(2)
contour(kxvalues,kyvalues,band2)
xlabel('k_x')
ylabel('k_y')
title('band 2')
axis('square')
