互注入半导体激光器的速率方程

clear;
i=sqrt(-1);
h=1e-12;
fs=800e9;%采样频率16GHz
Ts=1/fs;
t=40e-9:Ts:60e-9;
t1=40.15e9;%耦合延时时间
%初值
%G_10=15.48;%mw
x10=1;
y10=1;
n10=3;
%G_20=13.43 ;
x20=0.3;
y20=0.2;
n20=2;
Y1=zeros(3,length(t)); 
Y1(:,1)=[x10,y10,n10];
Y2=zeros(3,length(t)); 
Y2(:,1)=[x20,y20,n20];

%四阶龙格库塔
for j=1:length(t)-1
    %进行时延赋值
    id=j-round(t1*fs);%  t1*fs：延时点数
    if  id<=0
        x1_t1=0;
        y1_t1=0;
    else
        x1_t1=Y1(1,id);
        y1_t1=Y1(2,id); 
    end
    if  id<=0
        x2_t1=0;
        y2_t1=0;
    else
        x2_t1=Y2(1,id);
        y2_t1=Y2(2,id);
    end
    %解方程组
    k1=f(Y1(:,j),Y2(:,j),x1_t1,y1_t1,x2_t1,y2_t1);
    k2=f(Y1(:,j)+h*k1(1:3)/2,Y2(:,j)+h*k1(4:6)/2,x1_t1,y1_t1,x2_t1,y2_t1);
    k3=f(Y1(:,j)+h*k2(1:3)/2,Y2(:,j)+h*k2(4:6)/2,x1_t1,y1_t1,x2_t1,y2_t1);
    k4=f(Y1(:,j)+h*k3(1:3),Y2(:,j)+h*k3(4:6),x1_t1,y1_t1,x2_t1,y2_t1);
    Y1(:,j+1)=Y1(:,j)+h/6*(k1(1:3)+2*k2(1:3)+2*k3(1:3)+k4(1:3));
    Y2(:,j+1)=Y2(:,j)+h/6*(k1(4:6)+2*k2(4:6)+2*k3(4:6)+k4(4:6));
end
x1=Y1(1,:);
y1=Y1(2,:);
G=(x1+i*y1).*conj(x1+i*y1);
G_mean=mean(G);
G_1=G-G_mean;
x2=Y2(1,:);
y2=Y2(2,:);
H=(x2+i*y2).*conj(x2+i*y2);
H_mean=mean(H);
H_1=H-H_mean;

function F=f(Y1,Y2,x1_t1,y1_t1,x2_t1,y2_t1)
%参量
b=3;
J=1.222;
t1=40.15e9;
rc=5.36e11;
rs=5.96e9;
rn=7.53e9;
rp=1.91e11;
w02=0;%简化数值计算
%初始值
x1=Y1(1);
y1=Y1(2);
n1=Y1(3);
x2=Y2(1);
y2=Y2(2);
n2=Y2(3);
%参量
f1=-20e9;%失谐频率
z1=0.05;%耦合强度LD1注入LD2
z2=0.05;%耦合强度LD2注入LD1
%方程组
dx1_dt=0.5*(rc*rn*n1/(rs*J)-rp.*(x1^2+y1^2-1))*(x1+b*y1)+2*pi*f1*y1+z2*rc*(x2_t1*cos(w02*t1)-y2_t1*sin(w02*t1));
dy1_dt=0.5*(rc*rn*n1/(rs*J)-rp*(x1^2+y1^2-1))*(-b*x1+y1)-2*pi*f1*x1+z2*rc*(y2_t1*cos(w02*t1)+x2_t1*sin(w02*t1));
dn1_dt=-(rs+rn*(x1^2+y1^2))*n1-rs*J*(x1^2+y1^2-1)+rs*rp*J/rc*(x1^2+y1^2)*(x1^2+y1^2-1);

dx2_dt=0.5*(rc*rn*n2/(rs*J)-rp*(x2^2+y2^2-1))*(x2+b*y2)+z1*rc*(x1_t1*cos(w02*t1)-y1_t1*sin(w02*t1));
dy2_dt=0.5*(rc*rn*n2/(rs*J)-rp*(x2^2+y2^2-1))*(-b*x2+y2)+z1*rc*(y1_t1*cos(w02*t1)+x1_t1*sin(w02*t1));
dn2_dt=-(rs+rn*(x2^2+y2^2))*n2-rs*J*(x2^2+y2^2-1)+rs*rp*J/rc*(x2^2+y2^2)*(x2^2+y2^2-1);
F=[dx1_dt;dy1_dt;dn1_dt;dx2_dt;dy2_dt;dn2_dt];
end
matlab中这个程序是否存在物理逻辑问题，

我运行出来一直都是处于稳态（如下图），不知道程序那出了问题？