
% ================================================================
% John Kerl
% kerl at math dot arizona dot edu
% 2008-02-08
% ================================================================

% ----------------------------------------------------------------
% Euler:
% psi = psi + i * dt * (c Dpsi - V psi)
% t = t + dt

% ----------------------------------------------------------------
% Fourth-order Runge-Kutta ("RK4") approximation to
%   dy/dt = f(t,y)
% with
%   y(t_0) = y_0
% is
%   y_{n+1} = y_n + h/6(k1 + 2 k2 + 2 k3 + k4)
%   t_{n+1} = t_n + h
% where
%   k1 = f(t_n,       y_n)
%   k2 = f(t_n + h/2, y_n + k1 h/2)
%   k3 = f(t_n + h/2, y_n + k2 h/2)
%   k4 = f(t_n + h,   y_n + k3 h)

% ----------------------------------------------------------------
% SE:
%   i hbar dpsi/dt = -hbar^2 / 2m Dpsi + V psi
%   i dpsi/dt = -hbar / 2m Dpsi + V/hbar psi
%   dpsi/dt = i (hbar / 2m Dpsi - V/hbar psi)
%   dpsi/dt = f(psi)  -- run RK4's at each x.
%
%   y += y + dt/6(k1 + 2 k2 + 2 k3 + k4)
%   t += t + dt
% where
%   k1 = f(psi)
%   k2 = f(psi + k1 h/2)
%   k3 = f(psi + k2 h/2)
%   k4 = f(psi + k3 h)
% i.e.
%   k1 = i (hbar / 2m D(psi)          - V/hbar (psi))
%   k2 = i (hbar / 2m D(psi + k1*h/2) - V/hbar (psi + k1*h/2))
%   k3 = i (hbar / 2m D(psi + k2*h/2) - V/hbar (psi + k2*h/2))
%   k4 = i (hbar / 2m D(psi + k3*h)   - V/hbar (psi + k3*h))
% With hbar = 1 and mass --> c:
%   k1 = i (c D(psi)          - V (psi))
%   k2 = i (c D(psi + k1*h/2) - V (psi + k1*h/2))
%   k3 = i (c D(psi + k2*h/2) - V (psi + k2*h/2))
%   k4 = i (c D(psi + k3*h)   - V (psi + k3*h))

xlo = -10;
xhi =  10;
dx  =   0.5;

ylo = -10;
yhi =  10;
dy  =   0.5;

dt  =   0.01;
h   =   0.01;
c   =   1.0;

niter = 1000;

make_eps = 0

[x, y] = meshgrid(xlo:dx:xhi, ylo:dy:yhi);
psi0 = exp(-x.^2 + -y.^2);

V = 0.0 * x * y;
psi0 = psi0 / sqrt(sum(sum(conj(psi0).*psi0))); % normalize

psie = psi0;
psir = psi0;
t = 0;

for iter = 1:niter
	% Euler:
	psie = psie + i * dt * (c * del2(psie) - V.*psie);

	% RK4:
	k1 = i * (c * del2(psir)          - V .*  psir);
	k2 = i * (c * del2(psir + k1*h/2) - V .* (psir + k1*h/2));
	k3 = i * (c * del2(psir + k2*h/2) - V .* (psir + k2*h/2));
	k4 = i * (c * del2(psir + k3*h)   - V .* (psir + k3*h));
	psir = psir + dt/6*(k1 + 2*k2 + 2*k3 + k4);

	%disp(abs(psie));
	%disp(abs(psir));
	%disp(sqrt(sum(sum(conj(psie).*psie))));
	%disp(sqrt(sum(sum(conj(psir).*psie))));
	psie = psie / sqrt(sum(sum(conj(psie).*psie))); % normalize
	psir = psir / sqrt(sum(sum(conj(psir).*psir))); % normalize

	t = t + dt;
end

figure; surf(x,y,V); shading flat;
if make_eps
	set(gcf, 'PaperPositionMode', 'Manual');
	set(gcf, 'PaperPosition', [0 0 5.0 2.0])
	print -depsc se2dv
end

figure; surf(x,y,(abs(psi0))); shading flat;
if make_eps
	set(gcf, 'PaperPositionMode', 'Manual');
	set(gcf, 'PaperPosition', [0 0 5.0 2.0])
	print -depsc se2d0
end

figure; surf(x,y,(abs(psie))); shading flat;
if make_eps
	set(gcf, 'PaperPositionMode', 'Manual');
	set(gcf, 'PaperPosition', [0 0 5.0 2.0])
	print -depsc se2de
end

figure; surf(x,y,(abs(psir))); shading flat;
if make_eps
	set(gcf, 'PaperPositionMode', 'Manual');
	set(gcf, 'PaperPosition', [0 0 5.0 2.0])
	print -depsc se2dr
end

