%numerical solution to Schroedinger equation for
%rectangular potential well with infinite barrier energy

length = 1;	                     %length of well (nm)
npoints=200;					 %number of sample points
x=linspace(0,length,npoints+1)';	 %position of sample points
num_sol=199;			             %number of solutions sought
dx = x(2)-x(1);

v = zeros(size(x));              %potential (eV)

[energy,phi,H]=schrod(length,npoints,v,num_sol);		%call Schrod.

hbar=1.054571596;              %Planck's constant (x10^34 J s)
echarge=1.602176462;		   %electron charge (x10^19 C)
emass=9.10938188;              %electron mass (x10^31 kg)

% compute theoretical values for energy levels of particle in a box
n = num_sol:-1:1; n = fliplr(n); % flip depending on what schrod returns
Eth = n.^2*hbar^2*pi^2/(2*emass*length^2)/echarge;

% prepare some state
psi = zeros(size(x));
ind = x > 0.2 & x < 0.3; psi(ind) = x(ind) - 0.2;
ind = x >= 0.3 & x < 0.35; psi(ind) = 0.1;
ind = x >= 0.35 & x < 0.45; psi(ind) = 0.45 - x(ind);
psi = 1i*psi;
psi = psi/sqrt(psi'*psi); disp(['<psi|psi>=' num2str(psi'*psi)]);

c = phi'*psi; % decomposes psi into phi basis (phi = matrix)

xop = sparse(diag(x)); % x-operator
pop = diag(ones(npoints,1),1)-diag(ones(npoints,1),-1); % p-operator
pop = sparse(pop)/(2i*dx); %note p-op is actually k operator (k = p/hbar)

p2 = diag(ones(npoints,1),1)+diag(ones(npoints,1),-1)-2*diag(ones(npoints+1,1)); % p-operator
p2 = -sparse(p2)/(dx^2); %note p^2-op is actually k^2 operator

% print at most 10 first energy levels
if num_sol > 10, nval = 10; else nval = num_sol; end
for i=1:nval
    fprintf('eigenenergy (%d) = %f, theory = %f eV\n', i, energy(i), Eth(i));
end

figure(1); % display potential inside the box
plot(x,v,'b');xlabel('Distance (nm)'),ylabel('Potential energy, (eV)');

figure(2); % display comparison between theory and schrod.m
plot(n, Eth, '-', n, energy, 'o');
xlabel('eigenvalue #'); ylabel('energy (eV)'); legend('theory', 'schrod');

figure(3); % show probability density for psi
plot(x, abs(psi).^2, '.-'); ylabel('|\psi(x)|^2'); xlabel('Distance (nm)');

fprintf('checking decomposition norm: %g\n', c'*c);
fprintf('expectation energy: %g eV\n', c'*diag(energy)*c);
fprintf('expectation energy (another method): %g eV\n', psi(2:end-1)'*H*psi(2:end-1));
[v,i] = max(abs(c).^2);
fprintf('most probable outcome in energy measurement: energy(%g) = %g eV\n', i, energy(i));
x0 = psi'*xop*psi; p0 = psi'*pop*psi;
fprintf('expectation position: %g nm\n', x0);

Warning: For real symmetric problems, must have number of eigenvalues k < n.
Using eig(full(A)) instead. 
<psi|psi>=1
eigenenergy (1) = 0.376022, theory = 0.376030 eV
eigenenergy (2) = 1.503997, theory = 1.504120 eV
eigenenergy (3) = 3.383645, theory = 3.384271 eV
eigenenergy (4) = 6.014503, theory = 6.016482 eV
eigenenergy (5) = 9.395922, theory = 9.400753 eV
eigenenergy (6) = 13.527067, theory = 13.537084 eV
eigenenergy (7) = 18.406919, theory = 18.425476 eV
eigenenergy (8) = 24.034275, theory = 24.065928 eV
eigenenergy (9) = 30.407745, theory = 30.458440 eV
eigenenergy (10) = 37.525757, theory = 37.603012 eV
checking decomposition norm: 1
expectation energy: 6.52674 eV
expectation energy (another method): 6.52674 eV
most probable outcome in energy measurement: energy(2) = 1.504 eV
expectation position: 0.325 nm

Published with MATLAB® 7.10