% Script to simulate tunneling in 1-D
% Isaac Lenton
% August 2017

barrierWidth = 0.12;
barrierHeight = 2.1;

xmax = 10;
dt = 1.0e-2;                    % Time step size
N = 512;                        % Simulation spatial resolution
x = (-N/2:N/2-1)*2*xmax/N;      % Position axis
p = (-N/2:N/2-1)*2*pi/xmax;     % Momentum axis

% Calculate potential
potential = 0.25*x.^2 + barrierHeight*exp(-(x/barrierWidth).^2);

% Calculate initial wavefunction
wavefunction = sqrt(sqrt(1/(2*pi)))*exp(-0.25*(x+4).^2);

% Calculate the propogators
kprop = exp(-1i*dt*p.^2);
xprop = exp(-1i*dt*potential);

% Setup plot
figure;
[ax, hp, hw] = plotyy(x, potential, x, abs(wavefunction).^2);
axis(ax(1), [-xmax xmax 0 25]); set(ax(1), 'Ytick', 0:5:25);
axis(ax(2), [-xmax xmax 0 1]); set(ax(2), 'Ytick', 0:0.1:1);
legend('Potential', '|Wavefunction|^2');
xlabel('Position');

% Evolve the system and show the animation
while ishandle(hw)
    
    % Evolve the wavefunction
    kphi = fftshift(fft(fftshift(wavefunction))).*kprop;
    wavefunction = fftshift(ifft(fftshift(kphi))).*xprop;
    
    set(hw, 'ydata', abs(wavefunction).^2);
    pause(0.05);
    drawnow;
end