Preventing Eigenvalues from Crossing - Eigenshuffle Problem

I am calculating a bandstructure, which requires that I take eigenvalues along a series of points that span my space. I found the code for 'eigenshuffle' which is supposed to sort the eigenvalues in a consistent way (so that the slopes of each eigenvalue set are continuous). This works for me in areas where a crossing has occurred and the difference in energy (the y axis) between the two eigenvalue functions is relatively small. However, it does not seems to work in the long term. I have included a picture showing my problem, and the basic code for my eigenvalue problem.

The below code will run so long as the eigenshuffle function is loaded (which can be found on FileExchange). How can I get eigenshuffle to work on sorting these eigenvalues? Why isn't it working?

    sizee = 101;
    a=2.45;
    a1 = a* [sqrt(3)/2 1/2];
    a2 = a* [sqrt(3)/2 -1/2];
    yy = linspace(-0.9831,0.9831,sizee);
    xx = linspace(2.5541,-2.5541,sizee);
    [Kx Ky] = meshgrid(xx,yy);

    yv = [0 2.5541 2.5541 0 -2.5541 -2.5541];
    xv = [0.9831 0.4915 0.4916 -0.9831 -0.4916 -0.4915];

    IN = inpolygon(Kx,Ky,xv,yv);
    for m=1:sizee
        for n=1:sizee
            if IN(m,n)==1
                KxHEX(m,n) = Kx(m,n);
                KyHEX(m,n) = Ky(m,n);
            else
                KxHEX(m,n) = NaN;
                KyHEX(m,n) = NaN;
            end
        end
    end

    hbarSQM = 7.62;
    d1 = [a/sqrt(3) 0];
    d = sqrt(d1(1)^2 + d1(2)^2); 
    Es = -17.52; 
    Ep = -8.97; 
    VssS = (-1.40) * hbarSQM/(2*d^2);
    VspS = (1.84) * hbarSQM/(2*d^2);
    VppS = (3.24) * hbarSQM/(2*d^2);
    VppP = (-0.81) * hbarSQM/(2*d^2);

    Const = ones(sizee,sizee);
    Es = Es * Const;
    Ep = Ep * Const;

    %phase factors:
    R1 = [0 0];
    R2 = a/2 * [-sqrt(3) 1];
    R3 = -a/2 * [sqrt(3) 1];
    %dot products:
    i = sqrt(-1);
    kR1 = Kx*R1(1) + Ky*R1(2);
    kR2 = Kx*R2(1) + Ky*R2(2);
    kR3 = Kx*R3(1) + Ky*R3(2);
    g0 = 1 + exp(-i*kR1) + exp(-i*kR3);
    g1 = 1 - (1/2) * exp(-i*kR2) - (1/2) * exp(-i*kR3);
    g2 = sqrt(3)/2 * (exp(-i*kR2) - exp(-i*kR3));
    g3 = 1 + (1/4) * exp(-i*kR2) + (1/4) * exp(-i*kR3);
    g4 = (3/4) * (exp(-i*kR2) + exp(-i*kR3));
    g5 = sqrt(3)/4 * (-exp(-i*kR2) + exp(-i*kR3));

    H = zeros(8,8,sizee,sizee);

    for m=1:8
        if m==1 || m==5
            H(m,m,:,:) = Es;
        else
            H(m,m,:,:) = Ep;
        end
    end

    %off diagonals:
    H(5,1,:,:) = VssS * conj(g0); %row5, col1
    H(5,2,:,:) = VspS * conj(g1);
    H(5,3,:,:) = VspS * conj(g2);
    H(6,3,:,:) = (VppS + VppP) * conj(g5);
    H(6,2,:,:) = VppS*conj(g3) + VppP*conj(g5);
    H(7,3,:,:) = VppS*conj(g4) + VppP*conj(g3);
    H(8,4,:,:) = VppP * conj(g0);
    H(7,2,:,:) = H(6,3,:,:);
    H(7,1,:,:) = H(5,3,:,:);
    H(6,1,:,:) = H(5,2,:,:);

    for m=5:8
        for n=1:4
            H(n,m,:,:) = conj(H(m,n,:,:));
        end
    end

    vectors = zeros(8,8,sizee,sizee);
    energies = zeros(8,sizee,sizee);
    values = zeros(8,8,sizee,sizee);
    vectorsSquare = zeros(8,8,sizee,sizee);
    energiesSquare = zeros(sizee,sizee,8);

    for yy=1:sizee
        [vectArray valArray] = eigenshuffle(H(:,:,:,yy));
        energies(:,:,yy) = valArray;
    end

%NOW THE PLOTTED SYMMETRY AXES:
%go from K to gamma
gamma = round(sizee/2);

index = [0 gamma];

while index(1) < gamma
    index(1) = index(1) + 1;
    for m=1:8
        linA(m,index(1)) = energies(m,index(1),index(2));
    end
    linA(9,index(1)) = Kx(index(1),index(2));
    linA(10,index(1)) = Ky(index(1),index(2));
end

for m=1:8
    linB(m,(gamma-index(2)+1)) = energies(m,index(1),index(2));
end
linB(9,gamma-index(2)+1) = Kx(index(1),index(2));
linB(10,gamma-index(2)+1) = Ky(index(1),index(2));

%now from gamma to M
while index(2) > 1
    index(2) = index(2) - 1;
    for m=1:8
        linB(m,gamma-index(2)+1) = energies(m,index(1),index(2));
    end
    linB(9,gamma-index(2)+1) = Kx(index(1),index(2));
    linB(10,gamma-index(2)+1) = Ky(index(1),index(2));
end

o = 1;
for m=1:8
    linC(m,o) = energies(m,index(1),index(2));
end
linC(9,o) = Kx(index(1),index(2));
linC(10,o) = Ky(index(1),index(2));
o=o+1;

%now from M to K
while index(1) > round(gamma/2)
    index(1) = index(1) - 1;
    for m=1:8
        linC(m,o) = energies(m,index(1),index(2));
    end
    linC(9,o) = Kx(index(1),index(2));
    linC(10,o) = Ky(index(1),index(2));
    o=o+1;
end

colors = rand(3,8);

kaxis = [linA(10,:) linB(9,:) abs(linC(10,:))+linB(9,length(linB(9,:)))];
lin = [linA linB linC];
figure
for m=1:8
    plot(kaxis,lin(m,:), 'Color', colors(:,m)); hold on
end
title('Tight Binding Bands Along Symmetry Axes in Graphene')
xlabel('k values')
ylabel('Energy in eV')