diff --git a/PS2/q2c.asv b/PS2/q2c.asv
new file mode 100644
index 0000000..830619f
--- /dev/null
+++ b/PS2/q2c.asv
@@ -0,0 +1,76 @@
+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);
+
+% Find first N analytic eigenvalues
+E_col_analytic = (1/q) * (hbar^2 * pi^2 * n.*n) / (2*m*L^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('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('POSITION [m]'); ylabel('PROBABILITY DENSITY [1/m]');
+yticks([0.02 0.04 0.06 0.08 0.10 0.12]);
+legend('n=2');
+axis([0 1e-9 0 0.04]);
+
+%{
+% Plot numerical eigenvalues
+figure(2); clf; h = plot(n,E_col,'kx'); grid on;
+set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
+xlabel('EIGENVALUE NUMBER'); ylabel('ENERGY [eV]');
+axis([0 100 0 40]);
+
+% Add analytic eigenvalues to above plot
+
+hold on;
+plot(n,E_col_analytic,'k-');
+legend({'Numerical','Analytical'},'Location','northwest');
+%}
+
diff --git a/PS2/q2c.m b/PS2/q2c.m
index 87df436..9b69784 100644
--- a/PS2/q2c.m
+++ b/PS2/q2c.m
@@ -4,19 +4,20 @@ clear all;
 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  = 1e-10;                 %lattice constant
-x  = a * n;                 %x-coordinates
+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 = 0*x; %0 potential at all x
+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
@@ -33,23 +34,33 @@ 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_50 = eigenvectors(:,50);
+phi_2 = eigenvectors(:,2);
 
 % find the probability densities of position for 1st and 50th eigenvectors
 
 P_1 = phi_1 .* conj(phi_1);
-P_50 = phi_50 .* conj(phi_50);
+P_2 = phi_2 .* conj(phi_2);
 
 % Find first N analytic eigenvalues
 E_col_analytic = (1/q) * (hbar^2 * pi^2 * n.*n) / (2*m*L^2);
 
-% Plot the probability densities for 1st and 50th eigenvectors
+% Plot the probability densities for 1st and 2nd eigenvectors
 
-figure(1); clf; h = plot(x,P_1,'kx',x,P_50,'k-');
+figure(1); clf; h = plot(r,P_1,'k-');
 grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
-xlabel('POSITION [m]'); ylabel('PROBABILITY DENSITY [1/m]');
-legend('n=1','n=50');
+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]);
+
+%{
 % Plot numerical eigenvalues
 figure(2); clf; h = plot(n,E_col,'kx'); grid on;
 set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
@@ -61,5 +72,5 @@ axis([0 100 0 40]);
 hold on;
 plot(n,E_col_analytic,'k-');
 legend({'Numerical','Analytical'},'Location','northwest');
-
+%}
 
diff --git a/PS2/q2c_fig1_1s.png b/PS2/q2c_fig1_1s.png
new file mode 100644
index 0000000..c71b686
Binary files /dev/null and b/PS2/q2c_fig1_1s.png differ
diff --git a/PS2/q2c_fig2_2s.png b/PS2/q2c_fig2_2s.png
new file mode 100644
index 0000000..6584036
Binary files /dev/null and b/PS2/q2c_fig2_2s.png differ