混沌复洛伦兹系统反同步MATLAB实现

clc;
clear all;
close all;


m=200000;   

th1=10;th2=28;th3=8/3;

e=zeros（m5）;

x=zeros（m5）; x（1:）=[-1-2-3-4-5];
y=zeros（m5）; y（1:）=[12345]; 
                                         
dT=0.0002;

for k=1:m-1

    
%驱动系统
k1=th1*（-y（k1）+y（k3））; 
k2=th1*（-y（k1）-k1*dT/2+y（k3）+dT/2）; 
k3=th1*（-y（k1）-k2*dT/2+y（k3）+dT/2）;
k4=th1*（-y（k1）-k3*dT+y（k3）+dT）;
y（k+11）=y（k1）+dT*（k1+2*k2+2*k3+k4）/6;   

k1=th1*（-y（k2）+y（k4））; 
k2=th1*（-y（k2）-k1*dT/2+y（k4）+dT/2）; 
k3=th1*（-y（k2）-k2*dT/2+y（k4）+dT/2）;
k4=th1*（-y（k2）-k3*dT+y（k4）+dT）;
y（k+12）=y（k2）+dT*（k1+2*k2+2*k3+k4）/6;   
 
k1=th2*y（k1）-y（k1）*y（k5）-y（k3）; 
k2=th2*（y（k1）+dT/2）-（y（k1）+dT/2）*（y（k5）+dT/2）-（y（k3）+k1*dT/2）;
k3=th2*（y（k1）+dT/2）-（y（k1）+dT/2）*（y（k5）+dT/2）-（y（k3）+k2*dT/2）;
k4=th2*（y（k1）+dT）-（y（k1）+dT）*（y（k5）+dT）-（y（k3）+k3*dT）;
y（k+13）=y（k3）+dT*（k1+2*k2+2*k3+k4）/6;

k1=th2*y（k2）-y（k2）*y（k5）-y（k4）; 
k2=th2*（y（k2）+dT/2）-（y（k2）+dT/2）*（y（k5）+dT/2）-（y（k4）+k1*dT/2）;
k3=th2*（y（k2）+dT/2）-（y（k2）+dT/2）*（y（k5）+dT/2）-（y（k4）+k2*dT/2）;
k4=th2*（y（k2）+dT）-（y（k2）+dT）*（y（k5）+dT）-（y（k4）+k3*dT）;
y（k+14）=y（k4）+dT*（k1+2*k2+2*k3+k4）/6;  

k1=y（k1）*y（k3）+y（k2）*y（k4）-th3*y（k5）;
k2=（y（k1）+dT/2）*（y（k3）+dT/2）+（y（k2）+dT/2）*（y（k4）+dT/2）-th3*（y（k5）+k1*dT/2）;
k3=（y（k1）+dT/2）*（y（k3）+dT/2）+（y（k2）+dT/2）*（y（k4）+dT/2）-th3*（y（k5）+k2*dT/2）;
k4=（y（k1）+dT）*（y（k3）+dT）+（y（k2）+dT）*（y（k4）+dT）-th3*（y（k5）+k3*dT）;
y（k+15）=y（k5）+dT*（k1+2*k2+2*k3+k4）/6;    %+kk*（y（k3）-x（k3）-k1*dT/2）
 
%%响应系统
k1=th1*（-x（k1）+x（k3））-th1*（y（k3）+x（k3））; 
k2=th1*（-x（k1）-k1*dT/2+x（k3）+dT/2）-th1*（y（k3）+x（k3）+dT）; 
k3=th1*（-x（k1）-k2*dT/2+x（k3）+dT/2）-th1*（y（k3）+x（k3）+dT）;
k4=th1*（-x（k1）-k3*dT+x（k3）+dT）-th1*（y（k3）+x（k3）+2*dT）;
x（k+11）=x（k1）+dT*（k1+2*k2+2*k3+k4）/6;   

k1=th1*（-x（k2）+x（k4））-th1*（y（k4）+x（k4））; 
k2=th1*（-x（k2）-k1*dT/2+x（k4）+dT/2）-th1*（y（k4）+x（k4）+dT）; 
k3=th1*（-x（k2）-k2*dT/2+x（k4）+dT/2）-th1*（y（k4）+x（k4）+dT）;
k4=th1*（-x（k2）-k3*dT+x（k4）+dT）-th1*（y（k4）+x（k4）+2*dT）;
x（k+12）=x（k2）+dT*（k1+2*k2+2*k3+k4）/6;   
 
k1=th2*x（k1）-x（k1）*x（k5）-x（k3）+x（k1）*x（k5）+y（k1）*y（k5）-th2*（y（k1）+x（k1））; 
k2=th2*（x（k1）+dT/2）-（x（k1）+dT/2）*（x（k5）+dT/2）-（x（k3）+k1*dT/2）+（x（k1）+dT/2）*（x（k5）+dT/2）+（y（k1）+dT/2）*（y（k5）+dT/2）-th2*（y（k1）+x（k1）+dT）;
k3=th2*（x（k1）+dT/2）-（x（k1）+dT/2）*（x（k5）+dT/2）-（x（k3）+k2*dT/2）+（x（k1）+dT/2）*（x（k5）+dT/2）+（y（k1）+dT/2）*（y（k5）+dT/2）-th2*（y（k1）+x（k1）+dT）;
k4=th2*（x（k1）+dT）-（x（k1）+dT）*（x（k5）+dT）-（x（k3）+k3*dT）+（x（k1）+dT）*（x（k5）+dT）+（y（k1）+dT）*（y（k5）+dT）-th2*（y（k1）+x（k1）+2*dT）;
x（k+13）=x（k3）+dT*（k1+2*k2+2*k3+k4）/6;

%k1=th2*x（k2）-x（k2）*x（k5）-x（k4）+kk4*（y（k4）-x（k4））+x（k2）*x（k5）+y（k2）*y（k5）-th2*（y（k2）+x（k2））; 
k1=th2*x（k2）-x（k2）*x（k5）-x（k4）+x（k2）*x（k5）+y（k2）*y（k5）-th2*（y（k2）+x（k2））; 
k2=th2*（x（k2）+dT/2）-（x（k2）+dT/2）*（x（k5）+dT/2）-（x