% Plots of the two bound-state wavefunctions for the finite square well, 
% 0 <= x <= 1 and 2mV_0/hbar^2 = 16.
 
% x-axis parameters
xscale = 1;
xres = .001;
xleft = (-1:xres:0)*xscale;
xinside = (xres:xres:1)*xscale;
xright = (1+xres:xres:2)*xscale;
x = horzcat(xleft,xinside,xright);

% potential well 
V0 = 16; % in units of hbar^2/2m
V = horzcat(V0*ones(1,length(xleft)),zeros(1,length(xinside)),V0*ones(1,length(xright)));

% first bound state
k1 = 2.060; 
q1 = 3.429;
A = 1; 
B = A*exp(i*k1); 
C = A + B;
D = A*exp(q1)*(exp(i*k1)+1);
psi1left = C*exp(q1*xleft);
psi1inside = A*exp(i*k1*xinside) + B*exp(-i*k1*xinside);
psi1right = D*exp(-q1*xright);
psi1 = horzcat(psi1left,psi1inside,psi1right);

% second bound state
k2 = 3.791; 
q2 = 1.276;
A2 = 1; 
B2 = -A2*exp(i*k2); 
C2 = A2 + B2;
D2 = A*exp(q2)*(exp(i*k2)-1);
psi2left = C2*exp(q2*xleft);
psi2inside = A2*exp(i*k2*xinside) + B2*exp(-i*k2*xinside);
psi2right = D2*exp(-q2*xright);
psi2 = horzcat(psi2left,psi2inside,psi2right);

% plot wavefunction, offset by energy as fraction of well depth 
% (i.e. E/V0 = k^2/(2mV_0/hbar^2))
clf;
plot(x,k1^2+abs(psi1).^2,x,k2^2+abs(psi2).^2,x,V);
axis([xleft(1),xright(1000),0,V0*1.2]);
xlabel('position (arbitrary units)');
ylabel('probability density (arbitrary normalization)')
titlestring = sprintf('Bound-state probability densities for a finite square well\n (offset by energy as fraction of well depth)');
title(titlestring);


% E. Boyda / SMC Phys125 / Fall '17. 