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

% ----------------------------------------------------------------
% 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.2;
dt  =   0.01;
h   =   0.01;
c   =   1.0;

niter = 10000;

x = [xlo:dx:xhi];
V = 0.0 * x;

%% ----------------------------------------------------------------
mn = size(x);
n  = mn(2);
H  = disclap(n) + diag(V);
[eigenvectors,D] = eig(H);

%psi0 = exp(-x.^2);
psi0 = eigenvectors(:,1); % First eigenvector
%psi0 = eigenvectors(:,2); % Second eigenvector
%psi0 = 0.4*eigenvectors(:,1) + 0.6*eigenvectors(:,2);

psi0 = psi0 / sqrt(sum(conj(psi0).*psi0)); % normalize

%% ----------------------------------------------------------------
plot(x, abs(psi0),'b')

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

for iter = 1:niter

	% Euler:
	% psi = psi + i * dt * (del2(psi,dx) - V.*psi)
	% t = t + dt
	psie = psie + i * dt * (c * del2(psie,dx) - V'.*psie);

	% RK4 SE:
	% k1 = i * (c * del2(psi,dx)          - V *  psi)
	% k2 = i * (c * del2(psi + k1*h/2,dx) - V * (psi + k1*h/2))
	% k3 = i * (c * del2(psi + k2*h/2,dx) - V * (psi + k2*h/2))
	% k4 = i * (c * del2(psi + k3*h,dx)   - V * (psi + k3*h))
	% psi += dt/6(k1 + 2 k2 + 2 k3 + k4)

	k1 = i * (c * del2(psir,dx)          - V' .*  psir);
	k2 = i * (c * del2(psir + k1*h/2,dx) - V' .* (psir + k1*h/2));
	k3 = i * (c * del2(psir + k2*h/2,dx) - V' .* (psir + k2*h/2));
	k4 = i * (c * del2(psir + k3*h,dx)   - V' .* (psir + k3*h));
	psir = psir + dt/6*(k1 + 2*k2 + 2*k3 + k4);

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

	t = t + dt;
end
%disp(abs(psie));
%disp(abs(psir));
hold on;plot(x, abs(psie),'r')
hold on;plot(x, abs(psir),'g')


