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');