下面是一个模拟一维系统的 Matlab 代码。如何将其更改为具有两个子晶格的准一维二元链？

问：

在附图中，我附上了两个这样的晶格。（a） 用于一维晶格，（b） 用于准一维晶格。在（a）中，t1和t2是左右两侧的两个耦合器。在（b）中，我们有一些新的耦合，如k1、k2、k3，它们耦合了用蓝色和红色表示的两个子晶格（k3类似于t1/t2）。对于（a），附上有效的matlab代码。我的问题是如何更改此代码，以便我可以为 （b） 编写它。code_1 是主文件，int.m 是函数文件。这两个文件都应作为 .m 文件分别保存在同一文件夹中以运行。我的猜测是code_1，应该更改第 20 到 28 行。但是，问题是如何？非常感谢您的帮助。提前非常感谢。这些代码附在下面。

code_1.m：

% file: code1.m

% The code uses the function int.m
% parameter values should be  manually changed 

clear all
global M M1 tmax

% input values 
h=0.1;
k=-pi/2;mu=-0.05;
t1=exp(-h);t2=exp(h);
ene=2*cosh(mu+1i*k);
V1=exp(-1i*k-mu);V2=1/V1;

% Lattice size (left and right sides)
N1=50;N2=50;

N=N1+N2+1;
M=zeros(N,N);
M1=M;  
for n=1:N-1
if n > N1
M(n+1,n)=t1;M(n,n+1)=t2;
else
M(n+1,n)=t2;M(n,n+1)=t1;
end
end
M(1,1)=V1*1;
M(N,N)=V2*1;

sp=-2.*1i*cosh(mu+1i*k);
X=(1+(sp/2)^2)^.5;
if abs(sp/2-X)>1
X=-X;
disp(abs(sp/2-X))
end

nn=[0:1:N-1];
amp=4/(2+sp/X);
x0=amp.*exp(1i.*k.*nn+mu.*nn+h.*nn);
x0(N1+1:N)=x0(N1+1:N).*exp(-2*h*[0:1:N2]);

%tmax=150;
tmax=50;
dt=tmax/4999;
tspan=[0:dt:tmax];

Nt=max(size(tspan));
options=odeset('RelTol',1e-9,'AbsTol',1e-12);
[t0,y]=ode45('int',tspan,x0,options);

for j=1:max(size(t0))
potenza(j)=sum(abs(y(j,:).^2));
y1(j,:)=(y(j,:))/(potenza(j)).^.5;
end

figure(1)
plot(t0,log(potenza.^.5))

figure(2)
pcolor(abs(y))
pcolor([1:1:N],t0,abs(y1))
shading interp
colormap hot
colorbar

国际：

function der=int(t,x) % In this line, t = to and x = y in the main code

global M M1 tmax 

f=-1i*10*1;

if t<20
f=0;
end
if t>25
f=0;
end

der=-1i.*M*x-1i.*M1.*f*x;