clear all

L = 8;
amnu = 1;

r = [0:L^3-1];
nx = mod(r,L);
ny = mod((r-nx)/L,L);
nz = mod((r-ny*L-nx)/L^2,L);

numeigs = 33;

V = zeros(L^3,1);
c0 = 0;

%%%%%%%%%%%%%%%

w0 = 1.D0;
w1 = 1.D0;

% getting the coordinates for nearest neighbors in the x,y,z 
% directions (xp is +1 in the x-direction, xm is -1 in the 
% x-direction, and so on)

r_xp = nz*L^2          + ny*L          + mod(nx+1,L);
r_xm = nz*L^2          + ny*L          + mod(nx-1,L);

r_yp = nz*L^2          + mod(ny+1,L)*L + nx;
r_ym = nz*L^2          + mod(ny-1,L)*L + nx;

r_zp = mod(nz+1,L)*L^2 + ny*L          + nx;
r_zm = mod(nz-1,L)*L^2 + ny*L          + nx;

M = sparse([1:L^3],[1:L^3],3.0/amnu*w0);  
M = M + sparse([1:L^3],r_xp+1,-0.5/amnu*w1);  
M = M + sparse([1:L^3],r_xm+1,-0.5/amnu*w1);
M = M + sparse([1:L^3],r_yp+1,-0.5/amnu*w1); 
M = M + sparse([1:L^3],r_ym+1,-0.5/amnu*w1);  
M = M + sparse([1:L^3],r_zp+1,-0.5/amnu*w1); 
M = M + sparse([1:L^3],r_zm+1,-0.5/amnu*w1); 

M = M + sparse([1:L^3],[1:L^3],V);  

[vv,dd] = eigs(M,numeigs,'sa');

energies = diag(dd);

format long
energies

px = 2*pi/L;

[0; 1 - cos(px); 2*(1 - cos(px)); 3*(1 - cos(px)); 1 - cos(2*px)]



