2维势能中波函数模拟的代码（更新）
在 @PDE2718 同学的提醒下，改善了一下算法。
参考了如下文献，用了 Fourier 赝谱和（最简单的） Strang 方法分割时间演化算子。
Mechthild Thalhammer, Marco Caliari, Christof Neuhauser, High-order time-splitting Hermite and Fourier spectral methods, Journal of Computational Physics, Volume 228, Issue 3, 2009, Pages 822-832, ISSN 0021-9991, https://doi.org/10.1016/j.jcp.2008.10.008.
L = 20; % Interval Length
N = 501; % 设置 N*N 格点
x = linspace(-L/2,L/2-L/N,N)'; % Coordinate vector
dx = x(2) - x(1); % Coordinate step
NT = 50000; % No. of time steps
TF = 50; T = linspace(0,TF,NT); % Time vector
dT = T(2)-T(1); % Time step
I1 = speye(N,N);
% 动能部分
e1 = ones(N,1);
% % for even N
% K = [2*pi*(-N/2:-1)/L,2*pi*(0:N/2-1)/L];
% for odd N
K = [2*pi*(-N/2+1/2:-1)/L,2*pi*(0:N/2-1/2)/L];
Kx = kron(K,e1);
Ky = Kx';
UMmt = exp(-1i*dT*(Kx.^2+Ky.^2)/2); % exp(-i*(-0.5*Laplacian)*dT), 处于动量空间
% 势能部分
V=197;
a = 0.3;
U = -V*kron(...
spdiags(exp(-(x).^2/a^2),0,N,N),...
spdiags(exp(-(x).^2/a^2),0,N,N)...
); % 势能，可修改
G = 7;
l = 1.5;
Ui = -1i*G*(...
kron(...
spdiags(exp(-(x+0.5*L)/l)+exp(-(-x+0.5*L)/l),0,N,N),I1)+...
kron(...
I1,spdiags(exp(-(x+0.5*L)/l)+exp(-(-x+0.5*L)/l),0,N,N))...
); % 边界的虚势能，有时需要根据动量大小修改
Ut = reshape(diag(U+Ui),N,N); % 总势能
UVh = exp(-1i*Ut*dT/2); % exp(-i*U*dT/2), 处于位矢空间
% Gauss 波包
xf = kron(x,e1);
yf = kron(e1,x);
a = 2;% 控制波包大小
ko = 20;% (x向)动量中心
% 位置中心
xo = -4.2;
yo = 0;
% 沿 x 向传播的波包
psi0 = exp(-((xf-xo).^2+(yf-yo).^2)/(a^2)).*exp(1i*xf*ko);
psi0 = psi0/sqrt(psi0'*psi0); %归一化
psi2 = reshape(psi0,N,N);
% 时间演化，Strang 方法
% exp(-i*H*dT)=exp(-i*U*dT/2)*exp(-i*(-0.5*Laplacian)*dT)*exp(-i*U*dT/2)+O(dT^3)
for t = 1:NT
psi2 = full(UVh.*psi2);
psik = full(UMmt.*fftshift(fft2(psi2)));
psi2 = full(UVh.*ifft2(ifftshift(psik)));
rho1 = (abs(psi2).^2); % 密度
% Timenow = t*dT;
% 按需要加上别的命令，比如画图等
end
