在 @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