clear all

c(1) = 0.0;
c(2) = -0.4;
c(3) = -0.8;
c(4) = -1.2;
c(5) = -1.6;

ctarget = -3.0;

L = 20;
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);

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

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

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;

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

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

for nn = 1:size(c,2)
    [vv(:,nn),dd(:,nn)] = eigs(H0+c(nn)*H1,1,'sa');
end

[vsmall,notused1,notused2] = qr(vv,0);

for nc = 1:100
    ctarget(nc) = -0.1*nc;
    exact_energy(nc) = eigs(H0+ctarget(nc)*H1,1,'SA');

    for nn = 1:size(c,2)
        ECmatrix = vsmall(:,1:nn)'*(H0+ctarget(nc)*H1)*vsmall(:,1:nn);
        ECmatrix = (ECmatrix + ECmatrix')/2;
        EC(nc,nn) = eigs(ECmatrix,1,'SA');
    end
end    

hold on
for nEC = 1:size(c,2)
    plot(ctarget,EC(:,nEC),'x')
end
plot(ctarget,exact_energy,'o')
xlabel('c')
ylabel('E_0(c)')
xlim([-10.5 0.5])
ylim([-9 1])
hold off

