Who can help me in writing the following program: phase(i) = exp(1i*dot(k(:,index),R(:,i))) This gives you the phase you're going to multiply with the hopping parameter. You'd have to write a loop that runs from i = 1:3 in order to multiply each contribution to the hopping parameter such that the product of all 3 would get you the total hopping parameter. k and R here are matrices that contain generated wavevectors and the coordinates for the atoms, respectively. Es = -4.2; Ep = 1.715; orbitals = [Es Ep Ep Ep]; Vpp_sigma = 2.72; Vpp_pii = -0.72; Vss_sigma = -2.08; Vsp_sigma = 2.48; a = 3.86; a1 = (a/2) * [sqrt(3), -1, 0]; a2 = (a/2) * [sqrt(3), 1, 0]; coordinates = [4.45714 0 0; %B 2.22857 0 0.46152]; %A R1 = coordinates(2, :)-coordinates(1,:); R2 = R1 + a1; R3 = R1 + a2; R = [R1' R2' R3']; n = dot(R,[ 0 0 1])/norm(R); %direction cosine l = dot(R,[1 0 0])/norm(R); m = dot(R, [0 1 0])/norm(R); tss = Vss_sigma; txs = l*Vss_sigma; tys = m*Vss_sigma; tzs = n*Vss_sigma; txy = l*m*Vpp_sigma - l*m*Vpp_pii; txz = l*n*Vpp_sigma - l*n*Vpp_pii; tyz = m*n*Vpp_sigma - m*n*Vpp_pii; tzz = n^2 *Vpp_sigma+(1-n^2)*Vpp_pii; txx = l^2 *Vpp_sigma+(1-l^2)*Vpp_pii; tyy = m^2 *Vpp_sigma+(1-m^2)*Vpp_pii; These will be your matrix elements that you'll want. Since the Hamiltonian matrix takes a block-diagonal form, you can create an A block and B block and set the Hamiltonian equal to [A B; B A]. Once you run the loop to get all the hopping parameters correct you can simply plot the eigenvalues of the full Hamiltonian matrix across the desired path to obtain the full band structure.