clear; close all;
x0 = -17; t0 = 0; m = 1;
sig_x = 2; k0 = 4;
L = 4; V0 = 4;
mode = 1;
xmin = -50; xmax = 50; Nx = 2000;
int_xmin = x0-6*sig_x; int_xmax = x0+6*sig_x;
Nk = 500; kmax = k0+4/sig_x; kmin = max(kmax/Nk, k0-4/sig_x);
tmin = 0; tmax = 20; Nt = 200;
ax = [xmin, xmax, -0.5, 0.5];
psi = @(x) 1/(2*pi*sig_x^2)^0.25/sqrt(1 + 1i*t0/(2*m*sig_x^2))...
*exp(-(x-x0-k0*t0/m).^2/(2*sig_x)^2/(1 + 1i*t0/(2*m*sig_x^2)))...
.*exp(1i*k0*(x-k0*t0/(2*m)));
x = linspace(int_xmin, int_xmax, Nx); dx = (int_xmax-int_xmin)/(Nx-1);
k = linspace(kmin, kmax, Nk); dk = (kmax-kmin)/(Nk-1);
[g1, k1] = FFT(fftresize(psi(x),Nx*2), int_xmin, dx);
g = interp1(k1, g1, k);
figure; plot(k, abs(g), '.');
axis([kmin, kmax, 0, max(abs(g))*1.1]);
xlabel k; title 'FFT of \psi(x)';
E = k.^2/(2*m);
x = linspace(xmin, xmax, Nx);
t = linspace(tmin, tmax, Nt);
figure; set(gcf, 'Unit', 'Normalized', 'Position', [0.1, 0.1, 0.4, 0.3]);
for it = 1:Nt
psi1 = zeros(size(x));
if mode == 1
A = zeros(1,Nk); B = A;
for i = 1:Nk
A(i) = integral(@(x)psi(x).*FSB(x,E(i),L,m,V0,0), ...
int_xmin, int_xmax, 'RelTol', 1e-16);
B(i) = integral(@(x)psi(x).*FSB(x,E(i),L,m,V0,1), ...
int_xmin, int_xmax, 'RelTol', 1e-16);
psi1 = psi1 + dk*(A(i)*FSB(x,E(i),L,m,V0,0)...
+ B(i)*FSB(x,E(i),L,m,V0,1))*exp(-1i*E(i)*t(it));
end
else
for i = 1:Nk
C1 = zeros(1,Nk); C2 = C1;
C1(i) = integral(@(x)psi(x).*conj(FSB2(x,E(i),L,m,V0,1)), ...
xmin, xmax, 'RelTol', 1e-16);
C2(i) = integral(@(x)psi(x).*conj(FSB2(x,E(i),L,m,V0,-1)), ...
xmin, xmax, 'RelTol', 1e-16);
psi1 = psi1 + dk*(C1(i)*FSB2(x,E(i),L,m,V0,1) ...
+ C2(i)*FSB2(x,E(i),L,m,V0,-1))*exp(-1i*E(i)*t(it));
end
end
clf;
plot(x, real(psi1)); hold on; plot(x, imag(psi1));
plot([xmin, -L, -L, L, L, xmax], [0, 0, 0.4, 0.4, 0, 0], 'k');
axis(ax); xlabel x; title(['t = ' num2str(t(it), '%.2f')]);
saveas(gcf, [num2str(it) '.png']);
end