Nested For Loops and Plotting

Before I begin, I am an undergrad student with little experience in the field of solid-state physics, so bear with me here. I have posted before on the topic to no avail, but I have made significant changes and I think I'm considerably closer to my goal than I was prior. I'm trying to plot the energy levels of a 2D lattice, but I am having issues producing any meaningful results and I believe this is due to the nested for loops. The plot I am trying to replicate is present in the attached PDF in section 4, figure 5c. I started by mapping my 4D indices to 2D indices for ease of plotting. I also set all of the initial values to 1 for the same ease. I then initialized the Hamiltonian Operator to generate a tridiagonal matrix of m1*n1 and m2*n2 dimensions. The code currently produces the correct matrix, however I do not know how to factor in the alpha values which will be the x axis. Alpha is dependent on the magnetic field, B, which one can change. I know that B will have to be an array, so my first thought was to include B as 0:0.1:1 so I'd have enough values, but this runs into some size issues within the for loop. How do I proceed to produce valid results? Do a third for loop for B? Thanks.

% Initial Values
h_bar = 1;  mass = 1;   phi = 1;   B = 1; % 0:0.1:1;
a = 1;  t = 1;  imag = sqrt(-1);   M = 1;

% Matrix Dimensionality
m1 = 4;   m2 = 4;   n1 = 4;   n2 = 4;  N = 4;
i = (m1-1)*N+n1;  j = (m2-1)*N+n2;   % Equations given by Dr. Kunal Das

% Initialize Hamiltonian Operator
A = zeros(i,j);   % Preallocating Matrix
for i = 1:(m1*n1)   % 16x16 matrix has 16 elements in each row, 256 overall, etc
    for j = 1:(m2*n2)   % See above comment, 16 elements in each column now
       alpha = (B*a^2)/(phi);
       E = -2*cos((2*pi/N).*(M+alpha));   % Eq 24
       if i == j
           A(i,j) = E;
       elseif i == j + 1
           A(i,j) = t*exp(2*pi*imag*mass*alpha);  % Eq 28
       elseif i == j - 1
           A(i,j) = t*exp(-2*pi*imag*mass*alpha); % Eq 28
       elseif j == i + 1
           A(i,j) = t;     % Eq 29
       elseif j == i - 1
           A(i,j) = t;     % Eq 29
       else 
           A(i,j) = 0;   % Everything except tridiagonals are 0
       end  
    end
end

Es = eig(A);   Espec = Es';

% Plotting
figure(1); clf;
plot(alpha,Es,'.k');
xlabel('\alpha'); ylabel('E/|t|');  xlim([0,1]);  ylim([-4,4]);
title('Energy Spectra for Finite Sections of a Square Lattice: 4x4');

Answer 1

[Cris LaPierre]

I highly recommend you use Stephen Cobeldick's solution for creating a tridiagonal matrix (see here). It's much cleaner.

Avoid using i or j as loop counters, as these are the built in constants for imaginary numbers. This also means it is unnecessary to create your variable imag=sqrt(-1). Just use 1i instead.

exp(-2*pi*1i*mass*alpha)

You are setting limits that are preventing you from seeing the results (specifically your xlim).

Given your approach, I do think you need to loop through your values of B. However, with Stephen's suggestion, this reduces down to a single loop.

Running your code with those changes suggests you still have some issues with your calculations, but at least the plot appears. The trick I use is that MATLAB automatically treats each column of data as a series. I can then plot the entire Es matrix against alpha in a single plot command.

% Initial Values
mass = 1;   
phi = 1;   
B = 0:0.1:1;
a = 1;  
t = 1;  
M = 1;

% Matrix Dimensionality
m1 = 4;   
m2 = 4;   
n1 = 4;   
n2 = 4;  
N = 4;
K = (m1-1)*N+n1;  
L = (m2-1)*N+n2;   % Equations given by Dr. Kunal Das

% Initialize Hamiltonian Operator
alpha = (B*a^2)/(phi);
E = -2*cos((2*pi/N).*(M+alpha));   % Eq 24
b = t*exp(-2*pi*1i*mass*alpha); % Eq 28
c = t*exp(2*pi*1i*mass*alpha);  % Eq 28

for r = 1:length(B)
    A = diag(E(r)*ones(1,K)) + diag(b(r)*ones(1,K-1),1) + diag(c(r)*ones(1,K-1),-1);
    Es(:,r) = eig(A);
end

% Plotting
plot(alpha,Es,'.-k');
xlabel('\alpha'); ylabel('E/|t|');  ylim([-4,4]); xlim([0,1]);  
title('Energy Spectra for Finite Sections of a Square Lattice: 4x4');