This is what was on the matlab file:

%-----------------------------INSTRUCTIONS----------------------------------- % You must run se_solve_p_wells.m % Make sure that your current folder is the same folder where the above files are. % You can change the mass in the constants section on line 30 % To change the barrier height in the case of Double well, change the E % value in line 173. %----------------------------------------------------------------------------

%************************ METHOD & UNITS ************************************ % Matrix method used to solving time independent Schrodinger Equation % Energies values are in eV (electron volts) % Lengths are in nm (nanometers) % All other units are S.I. % Energy eignevalues displayed in Command Window % Energy spectrum and Well displayed in Figure Window %****************************************************************************

%*************************** POTENTIAL WELLS ******************************** % Different potential Wells for solving t-independent Schrodinger equation. %****************************************************************************

%---------------------------- CONSTANTS -------------------------------------

num = 801; % Number of data points hbar = 1.055e-34; % J.s e = 1.602e-19; % C me = 9.109e-31; % kg electron mass mp = 1.67252e-27; % kg proton mass mn = 1.67482e-27; % kg neutron mass eps0 = 8.854e-12; % F/m

m = me; % Mass of particle

Ese = 1.6e-19; % Energy scaling factor Lse = 1e-9; % Length scaling factor Cse = -hbar^2/(2*m) / (Lse^2*Ese); % Schrodinger Eq constant

%--------------------------------------------------------------------------

% Potential well parameters U = zeros(num,1); U_matrix = zeros(num-2);

% Potential well types disp(' '); disp(' '); disp(' POTENTIAL WELL TYPES - enter 1 or 2 or 3'); disp(' '); disp(' 1: Particle in a Box'); disp(' '); disp(' 2: Square Well'); disp(' '); disp(' 3: Sloping Particle in a Box'); % disp(' '); % disp(' 4: Harmonic Oscillator well'); % disp(' '); % disp(' 5: |x|'); % disp(' '); % disp(' 6: x^4'); % disp(' '); % disp(' 7: Double Well'); disp(' ');

wellType = input('Specify well type: 1, 2, 3 '); disp(' '); disp(' ');

% Potential Wells switch wellType %------------------------------Particle in a Box---------------------------

case 1 xMin = -0.1; xMax = +0.1; x1 = 0.05; % 1/2 well width U1 = 0; % Depth of well (eV) Ecut=20000; x = linspace(xMin,xMax, num); for cn = 1 : num if abs(x(cn)) = x1, U(cn) = Ecut; end end s = sprintf('Potential Well: Particle in a Box');

%------------------------------Square Well --------------------------------

case 2 xMin = -0.1; xMax = +0.1; U1 = 0; % Potential in the well x1 = 0.05; % 1/2 width of well Ecut=800; % Potential outside of the well x = linspace(xMin,xMax, num); for cn = 1 : num if abs(x(cn))=x1, U(cn) = Ecut; end; end s = sprintf('Potential Well: Square'); %-----------------------Sloping Particle in a Box--------------------------

case 3 xMin = -0.1; xMax = +0.1; U1 = 0; % Depth of LHS well U2 = 200; % Depth of RHS well x1 = 0.05; % 1/2 width of well Ecut=20000; x = linspace(xMin,xMax, num); intercept = (U1+U2)/2; slope = (U2-U1)/(2*x1); for cn = 1 : num if abs(x(cn))=x1, U(cn) = Ecut; end; end s = sprintf('Potential Well: Particle in a Sloping Box'); %----------------------------HARMONIC OSCILLATOR---------------------------

case 4 xMin = -0.4; xMax = +0.4; x1 = 0.2; % width U1 = 400; % well depth Ecut=2000; x = linspace(xMin,xMax, num); for cn = 1 : num if abs(x(cn))=x1/2, U(cn) = Ecut; end; end s = sprintf('Potential Well: HARMONIC OSCILLATOR'); %-----------------------------------MOD X----------------------------------

case 5 xMin = -0.15; xMax = +0.15; U1 = 400; % Depth of LHS well x1 = 0.1; % 1/2 width of well x = linspace(xMin,xMax, num); Ecut=400; slope = U1/x1; for cn = 1 : num if abs(x(cn))=x1, U(cn) = Ecut; end; end s = sprintf('Potential Well: |x|'); %----------------------------------X^4------------------------------------- case 6 xMin = -5.0; xMax = +5.0; x1 = 1.2; % 1/2 width of well Ecut= 2.07; x = linspace(xMin,xMax, num); for cn = 1 : num if abs(x(cn))=x1, U(cn) = Ecut; end end s = sprintf('Potential Well: x^4');

%------------------------------DOUBLE WELL--------------------------------- case 7 xMin = -5; xMax = +5; x1 =1.6; % 1/2 width of well x = linspace(xMin,xMax, num); A = 2.0; E = 2.0; % Energy Barrier Ecut=5.0; B = 2*sqrt(A*E); for cn = 1 : num if abs(x(cn))=x1, U(cn) = Ecut; end end s = sprintf('Potential Well: DOUBLE WELL'); end

% Potential energy matrix dx = (x(2)-x(1)); dx2 = dx^2;

for cn =1:(num-2) U_matrix(cn,cn) = U(cn+1); end

% Second Derivative Matrix off = ones(num-3,1); SD_matrix = (-2*eye(num-2) + diag(off,1) + diag(off,-1))/dx2;

% Kinetic Energy Matrix K_matrix = Cse * SD_matrix;

% Hamiltonian Matrix H_matrix = K_matrix + U_matrix;

% Eignevalues E_n and Eigenfunctions psi_N [e_funct , e_values] = eig(H_matrix);

% All Eigenvalues 1, 2 , ... n where E_N Ecut, flag = 1; end % if n = n + 1; end E(n-1) = []; n = n-2;

% Normalizing the wavefunction % Simpson's rule is used for integration for cn = 1 : n psi(:,cn) = [0; e_funct(:,cn); 0]; area = simpson1d((psi(:,cn) .* psi(:,cn))',xMin,xMax); psi(:,cn) = psi(:,cn)/sqrt(area); prob(:,cn) = psi(:,cn) .* psi(:,cn); if psi(5,cn)

% Display eigenvalues in Command Window disp(' '); disp('================================================================ '); disp(' '); fprintf('No. bound states found = %0.0g ',n); disp(' '); disp('Quantum State / Eigenvalues En (eV)');

for cn = 1 : n fprintf(' %0.0f ',cn); % state fprintf(' %0.5g ',E(cn)); % corresponding eigenvalue end disp(' ') disp(' ');

% Plot energy spectrum xs(1) = xMin; xs(2) = xMax;

figure(2); set(gcf,'Name','Energy Spectrum','NumberTitle','off') set(gcf,'color',[1 1 1 ]); set(gca,'fontSize',13);

plot(x,U,'b','LineWidth',2); xlabel('Position x (nm)','FontSize',14); ylabel('Energy U, E_n (eV)','FontSize',14); h_title = title(s); set(h_title,'FontSize',12); hold on

%cnmax = length(E); cnmax = 5

for cn = 1 : cnmax ys(1) = E(cn); ys(2) = ys(1); plot(xs,ys,'r','LineWidth',2); end axis([xMin-eps xMax min(U) E(5)+200]);

% Plots first 5 wavefunctions & probability density functions if n

figure(11) clf set(gcf,'NumberTitle','off'); set(gcf,'Name','Eigenvectors & Prob. densities'); set(gcf,'Color',[1 1 1]);

for cn = 1:nMax subplot(nMax,2,2*cn-1); y1 = psi(:,cn) ./ (max(psi(:,cn)-min(psi(:,cn)))); y2 = 1 + 2*U ./ (max(U)-min(U)); plot(x,y1,'lineWidth',2) hold on plot(x,y1(1)*ones(size(x)),'k--','linewidth',2) plot(x,y2-2,'r','lineWidth',1) axis off title_m = ['\psi n = ', num2str(cn)] ; title(title_m,'Fontsize',10); subplot(nMax,2,2*cn); y1 = prob(:,cn) ./ max(prob(:,cn)); y2 = 1 + 2 * U ./ (max(U) - min(U)); plot(x,y1,'lineWidth',2) hold on plot(x,y1(1)*ones(size(x)),'k--','linewidth',2) plot(x,y2-2,'r','lineWidth',1) title_m = ['\psi^2 n = ', num2str(cn)] ; title(title_m,'Fontsize',10); axis off end function integral = simpson1d(f,a,b)

% [1D] integration - Simpson's 1/3 rule % f function a = lower bound b = upper bound % Must have odd number of data points % Simpson's coefficients 1 4 2 4 ... 2 4 1

numS = length(f); % number of data points

sc = 2*ones(numS,1); sc(2:2:numS-1) = 4; sc(1) = 1; sc(numS) = 1;

h = (b-a)/(numS-1);

integral = (h/3) * f * sc; end

In class we have discussed in detail the particle (electron) in a box as an approximate model for an electron in a one-dimensional metal nanoparticle. In this exercise you will compare the results for this model with those for the square well model which is a more accurate representation of the system in that it includes the finite energy cost for removing the electron from the particle (the ionization potential or work function). In addition, you will consider the particle in a box model with a sloping potential in the box, which is representative of the case where an electric field is applied to the particle. Use the attached Matlab code, which numerically calculates the energy eigenvalues and eigen- functions, to answer the following questions. Put your answers in a Word or PDF document, typed or a scanned handwritten page. 1) For the particle in a box model, calculate the four lowest energy levels for the confined electron and show that they are consistent with the energy eigenvalue equation we derived in class. Deter- mine the length of the box in using the ground state (n 1) energy. In class we have discussed in detail the particle (electron) in a box as an approximate model for an electron in a one-dimensional metal nanoparticle. In this exercise you will compare the results for this model with those for the square well model which is a more accurate representation of the system in that it includes the finite energy cost for removing the electron from the particle (the ionization potential or work function). In addition, you will consider the particle in a box model with a sloping potential in the box, which is representative of the case where an electric field is applied to the particle. Use the attached Matlab code, which numerically calculates the energy eigenvalues and eigen- functions, to answer the following questions. Put your answers in a Word or PDF document, typed or a scanned handwritten page. 1) For the particle in a box model, calculate the four lowest energy levels for the confined electron and show that they are consistent with the energy eigenvalue equation we derived in class. Deter- mine the length of the box in using the ground state (n 1) energy