A question about the matlab error "Subscript indices must either be real positive integers or logicals"

Question

I wrote a simple matlab code to evaluate the ground state energy of a set of tight-binding Hamiltonians with local magnetic impurities. I always encounter the subscript error until I add a new line, "clear H1u H1d eigu eigd eig", before the end of the for loop. (See the attached code below).

I figure this out by trial and error but I have no idea why it works. Can anybody help me to clarify the problem?

Many thanks.

ps. I also notice that if I use, say "eigs(H1u, 3*N, 'SA')", in evaluating the lowest 3N eigenvalues of the corresponding matrix, the same error will not appear and the code works just fine. I also don't know why.

------------The Matlab code ----------------

clear all
%% 1D spinless tight-binding diamond chain
N=10; %% number of unit cells
h=zeros(3*N, 3*N); %%% the size of the Hamiltonian
h(1, 3*N)=-1; h(1, 3*N-1)=-1; h(1, 2)=-1; h(1, 3)=-1; %% specify periodic bc
for n=2:N
    h(3*n-2, 3*n-1)=-1; 
    h(3*n-2, 3*n)=-1;
    h(3*n-4, 3*n-2)=-1;
    h(3*n-3, 3*n-2)=-1;
end
h=h+h';  %% the complete tight-binding Hamiltonian
%------------------------------
J=0.2;  %% exchange coupling
L=5;  %% the maximum impurity separation

% Calculate the ground state energy of the half-filled band using eig
Epp=zeros(L, 1);

for n=1:L
H1u=h;  %% spin up Hamiltonian
H1d=h;  %% spin down Hamiltonian 
H1u(1, 1)=0.5*J;
H1u(3*n+1, 3*n+1)=0.5*J;
H1d(1, 1)=-0.5*J;
H1d(3*n+1, 3*n+1)=-0.5*J;
eigu=eig(H1u);
eigd=eig(H1d);
eig=[eigu; eigd];
eig=sort(eig); %% sort the single-particle energy levels
Epp(n)=sum(eig(1:3*N)); %%% the ground state energy of the half-filled band
clear H1u H1d eigu eigd eig 
end

Erkky=Epp;

plot(1:L, Erkky,'s-k');
xlabel('Impurity separation','FontSize',12);
ylabel('\Delta E','FontSize',12)

Answer 1

eig is the name of a Matlab function, but it's also the name of one of your variables.

If you redefine eig to be an array, as you do in the line eig = [eigu; eigd], then the line eigu = eig (H1u) no longer works.

Answer 2

Change the following codes

eig=[eigu; eigd];
eig=sort(eig); %% sort the single-particle energy levels
Epp(n)=sum(eig(1:3*N)); %%% the ground state energy of the half-filled band

as

eig_s=[eigu; eigd];
eig_s=sort(eig_s); %% sort the single-particle energy levels
Epp(n)=sum(eig_s(1:3*N)); %%% the ground state energy of the half-filled band