clear all;
%physical constants in MKS units

hbar = 1.054e-34;
q = 1.602e-19;
m = 9.110e-31;
epsilon_0 = 8.854e-12;

%generate lattice

N  = 100;                   %number of lattice points
n  = [1:N];                 %lattice points
a  = 0.1e-10;                 %lattice constant
r  = a * n;                 %x-coordinates
t0 = (hbar^2)/(2*m*a^2)/q;  %encapsulating factor
L  = a * (N+1);             %total length of consideration

%set up Hamiltonian matrix 

U = -q^2./(4*pi*epsilon_0.*r) * (1/q); %potential at r in [eV]
main_diag  = diag(2*t0*ones(1,N)+U,0); %create main diagonal matrix
lower_diag = diag(-t0*ones(1,N-1),-1); %create lower diagonal matrix
upper_diag = diag(-t0*ones(1,N-1),+1); %create upper diagonal matrix

H = main_diag + lower_diag + upper_diag; %sum to get Hamiltonian matrix

[eigenvectors,E_diag] = eig(H); %"eigenvectors" is a matrix wherein each column is an eigenvector
                                %"E_diag" is a diagonal matrix where the
                                %corresponding eigenvalues are on the
                                %diagonal.
                                
E_col = diag(E_diag); %folds E_diag into a column vector of eigenvalues

% return eigenvectors for the 1st and 50th eigenvalues

phi_1  = eigenvectors(:,1);
phi_2 = eigenvectors(:,2);

% find the probability densities of position for 1st and 50th eigenvectors

P_1 = phi_1 .* conj(phi_1);
P_2 = phi_2 .* conj(phi_2);

% Plot the probability densities for 1st and 2nd eigenvectors

figure(1); clf; h = plot(r,P_1,'k-');
grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('RADIAL POSITION [m]'); ylabel('PROBABILITY DENSITY [1/m]');
yticks([0.02 0.04 0.06 0.08 0.10 0.12]);
legend('n=1');
axis([0 1e-9 0 0.12]);

figure(2); clf; h = plot(r,P_2,'k-');
grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('RADIAL POSITION [m]'); ylabel('PROBABILITY DENSITY [1/m]');
yticks([0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04]);
legend('n=2');
axis([0 1e-9 0 0.04]);