模擬Feshbach的視頻中用到的部分代碼：

% 設(shè)置格點化的空間坐標(biāo)

L = 50; % Interval Length

N = 2000; % No of points

x = linspace(0.5,0.5+L,N)'; % Coordinate vector

dx = x(2) - x(1); % Coordinate step

%***************************************************************

%設(shè)置 Hamiltonian

e = ones(N,1);

e2 = [1;0];% open channel對應(yīng)的（二能級）態(tài)矢量

e3 = [0;1];% close channel對應(yīng)的（二能級）態(tài)矢量

e4 = [1;1];

V = 10;% 控制相互作用的強(qiáng)度

Lap = spdiags([e -2*e e],[-1 0 1],N,N)/dx^2; % 1維Laplacian

Uclose = 1e3./(x.^12)-1e3./(x.^6)+ 107.965; % close channel 的勢能

Uopen = (1e3./(x.^12)-1.5e3./(x.^6));% open channel 的勢能

C = V*(pi/2-atan(x-1.5));% 在左側(cè)加上channel間的耦合

G =0.5;% 控制邊界上的衰減

U1 = -(1i*G).*(atan(0.5*(x-(L)))+pi/2).*(1./(1+(3e-2)*(L-x).^2));

% 在右側(cè)緩緩加上虛勢能，讓反射的低能量的波包也可以被移出系統(tǒng)

% plot(x,Uopen,x,Uclose)

% axis([0.5 10 -150 105])

% Hamiltonian

H1 = -(1/2)*kron(Lap,speye(2)) +...? ? ? % 動能

? ? kron(spdiags(Uclose,0,N,N),spdiags(e3,0,2,2))+... %close channel potential

? ? kron(spdiags(Uopen,0,N,N),spdiags(e2,0,2,2))+... %open channel potential

? ? kron(spdiags(C,0,N,N),spdiags([e4 e4],[-1 1],2,2))+... % coupling between channels

? ? kron(spdiags(U1,0,N,N),speye(2)); % 虛勢能

% Parameters for computing psi(x,T) at different times 0 < T < TF

NT = 1000; % No. of time steps

TF = 80; T = linspace(0,TF,NT); % Time vector

dT = T(2)-T(1); % Time step

hbar = 1;

% Time displacement operator E=exp(-iHdT/hbar)

E = expm(-1i*full(H1)*dT/hbar); % time desplacement operator

% ********************

%設(shè)置初態(tài)

ko1 = -0.5; % Peak momentum

a = 4; % package width

xo=15;% peak center

% 離散化的空間坐標(biāo)對最大動量有限制，不然會出現(xiàn)Bloch振蕩這類現(xiàn)象

if ko1>(pi*N/L)

? ? error('too large ko1');

end

S1 = [1;0];

S1 = S1/sqrt(S1'*S1);% 設(shè)置初態(tài)channel

% 設(shè)置空間波函數(shù)為 平面波*Gauss輪廓

psi1 = exp(-(x-xo).^2/(2*a.^2)).*exp(1i*x*ko1);

psi1 = psi1/sqrt(psi1'*psi1*dx); % normalize the psi(x,0)

% 把兩個通道的自由度和空間自由度的態(tài)矢量張量積起來

psi1 = kron(psi1,S1);