% Initialize Parameters
N = 500; % Number of grid points
L = 100; % Width of axis
x = linspace(-L,L,N); % Grid points
h = (2*L) / N; % Grid point increments
h_bar = 1; mass = 1; w = 1; % Natural Units
tau = 1;
% Potential Barrier
V0 = 0.1;
V = zeros(1,N);
a = 5;
bvals = linspace(-8,8,N/4);
num_b = length(bvals);
for bidx = 1 : num_b % N divided by 4 for speed purposes
b = bvals(bidx);
for i = 1:N
if x(i) <= a+b && x(i) >= -a+b
V(i) = V0;
end
end
% Set up the Hamiltonian Operator Matrix
ham = zeros(N); % Set all elements to zero
coeff = -h_bar^2/(2*mass*h^2);
for i = 2:(N-1)
ham(i,i-1) = coeff;
ham(i,i) = -2*coeff + V(i);
ham(i,i+1) = coeff;
end
% First and last rows for periodic boundary conditions
ham(1,N) = coeff; ham(1,1) = -2*coeff; ham(1,2) = coeff;
ham(N,N-1) = coeff; ham(N,N) = -2*coeff; ham(N,1) = coeff;
% Compute the Crank-Nicolson matrix
dCN = ( inv(eye(N) + .5*j*tau/h_bar*ham) * (eye(N) - .5*j*tau/h_bar*ham));
% Initialize the wavefunction
x0 = 0; % Location of the center of the wavepacket
velocity = 0.5;
k0 = mass*velocity/h_bar; % Wavenumber
sigma0 = 5; % Width of wavefunction
Norm_psi = 1/(sqrt(sigma0*sqrt(pi))); % Normalization
psi = Norm_psi * cos(k0*x') .* exp(-(x'-x0).^2/(2*sigma0^2));
% Time Evolution of Psi
max_iter = L/(2*velocity*tau);
for iter = 1:max_iter
psi = dCN*psi;
end
% Transmission
xR = x;
xL = x;
for i = 1:length(x)
if x(i)>= a+b
xL(i) = 0;
end
if x(i)<= -a+b
xR(i) = 0;
end
end
pFinal = psi.*conj(psi);
tR = trapz(xR,pFinal);
tL = trapz(xL,pFinal);
transmission(bidx) = abs(tR - tL)/2;
end
% Plot probability versus position at various times
figure(1)
plot(x,pFinal);
xlabel('x'); ylabel('P(x,t)');
title('Probability density at various times');
r = rectangle('Position',[(-a+b) 0 2*(a) max(max(pFinal))]);
r.LineStyle = '-.';
plot(bvals, transmission); xlabel('b'); ylabel('transmission')