%This is the program for Generalized synchronization of fractional order Hyperchaotic chen system, a simple but fast
%method in time domain.
clc
clear
clf
a=35;b=3;c=12;d=7;r=0.5; %parameter a b c
q1=0.95;q2=0.95;q3=0.95;q4=0.95; %the fractional order
%q1=0.96;q2=0.96;q3=0.96;q4=0.96;
%q1=0.97;q2=0.97;q3=0.97;q4=0.97;
%q1=0.98;q2=0.98;q3=0.98;q4=0.98;
%q1=0.99;q2=0.99;q3=0.99;q4=0.99;
%q1=1;q2=1;q3=1;q4=1;
%q1=0.94;q2=0.94;q3=0.94;q4=0.94;


x10=2;x20=3;x30=1.5;x40=1.2;% initial values 
y10=2.1;y20=3.2;y30=1.2;y40=1.2;
z10=2;z20=3;z30=1;z40=1;

h=0.01;N=1000;
t1=0;t(N+1)=[0];

x1(N+1)=[0];x2(N+1)=[0];x3(N+1)=[0];x4(N+1)=[0];
xx1(N+1)=[0];xx2(N+1)=[0];xx3(N+1)=[0];xx4(N+1)=[0];
y1(N+1)=[0];y2(N+1)=[0];y3(N+1)=[0];y4(N+1)=[0];
yy1(N+1)=[0];yy2(N+1)=[0];yy3(N+1)=[0];yy4(N+1)=[0];
z1(N+1)=[0];z2(N+1)=[0];z3(N+1)=[0];z4(N+1)=[0];
zz1(N+1)=[0];zz2(N+1)=[0];zz3(N+1)=[0];zz4(N+1)=[0];

M1=0;M2=0;M3=0;M4=0;
M5=0;M6=0;M7=0;M8=0;
M9=0;M10=0;M11=0;M12=0;

N1=0;N2=0;N3=0;N4=0;
N5=0;N6=0;N7=0;N8=0;
N9=0;N10=0;N11=0;N12=0;

xx1(1)=x10+h^q1*(a*(x20-x10)+x40)/(q1*gamma(q1));
xx2(1)=x20+h^q2*((d-x30)*x10+c*x20)/(q2*gamma(q2));
xx3(1)=x30+h^q3*(x10*x20-x30*b)/(q3*gamma(q3));
xx4(1)=x40+h^q4*(x20*x30+r*x40)/(q4*gamma(q4));

yy1(1)=y10+h^q1*(-35*y10+17.5*y20-0.5*y40)/(q1*gamma(q1));
yy2(1)=y20+h^q2*(-64.062*y10+31.9723*y20+0.9706*(x40-y40)+x10*(x30-39.031)+19.9723*x20)/(q2*gamma(q2));
yy3(1)=y30+h^q3*(-y30-2*x30+x10*x20)/(q3*gamma(q3));
yy4(1)=y40+h^q4*(1.8303*y10-0.9706*y20-1.9723*y40+0.9152*x10+x20*(x30-0.9706)+2.4723*x40)/(q4*gamma(q4));

zz1(1)=z10+h^q1*(-35*z10+17.5*z20-0.5*z40)/(q1*gamma(q1));
zz2(1)=z20+h^q2*(-64.062*z10+31.9723*z20+0.9706*(x40-z40)+x10*(x30-39.031)+19.9723*x20)/(q2*gamma(q2));
zz3(1)=z30+h^q3*(-z30-2*x30+x10*x20)/(q3*gamma(q3));
zz4(1)=z40+h^q4*(1.8303*z10-0.9706*z20-1.9723*z40+0.9152*x10+x20*(x30-0.9706)+2.4723*x40)/(q4*gamma(q4));

x1(1)=x10+h^q1*(a*(xx2(1)-xx1(1))+xx4(1)+q1*(a*(x20-x10)+x40))/gamma(q1+2);
x2(1)=x20+h^q2*((d-xx3(1))*xx1(1)+c*xx2(1)+q2*((d-x30)*x10+c*x20))/gamma(q2+2);
x3(1)=x30+h^q3*(xx1(1)*xx2(1)-xx3(1)*b+q3*(x10*x20-x30*b))/gamma(q3+2);
x4(1)=x40+h^q4*(xx2(1)*xx2(1)+xx4(1)*r+q4*(x20*x30+x40*r))/gamma(q4+2);

y1(1)=y10+h^q1*(-35*yy1(1)+17.5*yy2(1)-0.5*yy4(1)+q1*(-35*y10+17.5*y20-0.5*y40))/gamma(q1+2);
y2(1)=y20+h^q2*(-64.062*yy1(1)+31.9723*yy2(1)+0.9706*(xx4(1)-yy4(1))+xx1(1)*(xx3(1)-39.031)+19.9723*xx2(1)+q2*(-64.062*y10+31.9723*y20+0.9706*(x40-y40)+x10*(x30-39.031)+19.9723*x20))/gamma(q2+2);
y3(1)=y30+h^q3*(-yy3(1)-2*xx3(1)+xx1(1)*xx2(1)+q3*(-y30-2*x30+x10*x20))/gamma(q3+2);
y4(1)=y40+h^q4*(1.8303*yy1(1)-0.9706*yy2(1)-1.9723*yy4(1)+0.9152*xx1(1)+xx2(1)*(xx3(1)-0.9706)+2.4723*xx4(1)+q4*(1.8303*y10-0.9706*y20-1.9723*y40+0.9152*x10+x20*(x30-0.9706)+2.4723*x40))/gamma(q4+2);

z1(1)=z10+h^q1*(-35*zz1(1)+17.5*zz2(1)-0.5*zz4(1)+q1*(-35*z10+17.5*z20-0.5*z40))/gamma(q1+2);
z2(1)=z20+h^q2*(-64.062*zz1(1)+31.9723*zz2(1)+0.9706*(xx4(1)-zz4(1))+xx1(1)*(xx3(1)-39.031)+19.9723*xx2(1)+q2*(-64.062*z10+31.9723*z20+0.9706*(x40-z40)+x10*(x30-39.031)+19.9723*x20))/gamma(q2+2);
z3(1)=z30+h^q3*(-zz3(1)-2*xx3(1)+xx1(1)*xx2(1)+q3*(-z30-2*x30+x10*x20))/gamma(q3+2);
z4(1)=z40+h^q4*(1.8303*zz1(1)-0.9706*zz2(1)-1.9723*zz4(1)+0.9152*xx1(1)+xx2(1)*(xx3(1)-0.9706)+2.4723*xx4(1)+q4*(1.8303*z10-0.9706*z20-1.9723*z40+0.9152*x10+x20*(x30-0.9706)+2.4723*x40))/gamma(q4+2);


for n=1:N
    t(n+1)=t(n)+h;
    M1=(n^(q1+1)-(n-q1)*(n+1)^q1)*(a*(x20-x10)+x40);
    M2=(n^(q2+1)-(n-q2)*(n+1)^q2)*((d-x30)*x10+c*x20);
    M3=(n^(q3+1)-(n-q3)*(n+1)^q3)*(x10*x20-x30*b);
    M4=(n^(q4+1)-(n-q4)*(n+1)^q4)*(x20*x30+x40*r);
    
    M5=(n^(q1+1)-(n-q1)*(n+1)^q1)*(-35*y10+17.5*y20-0.5*y40);
    M6=(n^(q2+1)-(n-q2)*(n+1)^q2)*(-64.062*y10+31.9723*y20+0.9706*(x40-y40)+x10*(x30-39.031)+19.9723*x20);
    M7=(n^(q3+1)-(n-q3)*(n+1)^q3)*(-y30-2*x30+x10*x20);
    M8=(n^(q4+1)-(n-q4)*(n+1)^q4)*(1.8303*y10-0.9706*y20-1.9723*y40+0.9152*x10+x20*(x30-0.9706)+2.4723*x40);
    
    M9=(n^(q1+1)-(n-q1)*(n+1)^q1)*(-35*z10+17.5*z20-0.5*z40);
    M10=(n^(q2+1)-(n-q2)*(n+1)^q2)*(-64.062*z10+31.9723*z20+0.9706*(x40-z40)+x10*(x30-39.031)+19.9723*x20);
    M11=(n^(q3+1)-(n-q3)*(n+1)^q3)*(-z30-2*x30+x10*x20);
    M12=(n^(q4+1)-(n-q4)*(n+1)^q4)*(1.8303*z10-0.9706*z20-1.9723*z40+0.9152*x10+x20*(x30-0.9706)+2.4723*x40);
    
    N1=((n+1)^q1-n^q1)*(a*(x20-x10)+x40);
    N2=((n+1)^q2-n^q2)*((d-x30)*x10+c*x20);
    N3=((n+1)^q3-n^q3)*(x10*x20-x30*b);
    N4=((n+1)^q4-n^q4)*(x20*x30+x40*r);
    
    N5=((n+1)^q1-n^q1)*(-35*y10+17.5*y20-0.5*y40);
    N6=((n+1)^q2-n^q2)*(-64.062*y10+31.9723*y20+0.9706*(x40-y40)+x10*(x30-39.031)+19.9723*x20);
    N7=((n+1)^q3-n^q3)*(-y30-2*x30+x10*x20);
    N8=((n+1)^q4-n^q4)*(1.8303*y10-0.9706*y20-1.9723*y40+0.9152*x10+x20*(x30-0.9706)+2.4723*x40);
    
    N9=((n+1)^q1-n^q1)*(-35*z10+17.5*z20-0.5*z40);
    N10=((n+1)^q2-n^q2)*(-64.062*z10+31.9723*z20+0.9706*(x40-z40)+x10*(x30-39.031)+19.9723*x20);
    N11=((n+1)^q3-n^q3)*(-z30-2*x30+x10*x20);
    N12=((n+1)^q4-n^q4)*(1.8303*z10-0.9706*z20-1.9723*z40+0.9152*x10+x20*(x30-0.9706)+2.4723*x40);
    
    for j=1:n
        M1=M1+((n-j+2)^(q1+1)+(n-j)^(q1+1)-2*(n-j+1)^(q1+1))*(a*(x2(j)-x1(j))+x4(j));
        M2=M2+((n-j+2)^(q2+1)+(n-j)^(q2+1)-2*(n-j+1)^(q2+1))*((d-x3(j))*x1(j)+c*x2(j));
        M3=M3+((n-j+2)^(q3+1)+(n-j)^(q3+1)-2*(n-j+1)^(q3+1))*(x1(j)*x2(j)-x3(j)*b);
        M4=M4+((n-j+2)^(q4+1)+(n-j)^(q4+1)-2*(n-j+1)^(q4+1))*(x3(j)*x2(j)+x4(j)*r);
        
        M5=M5+((n-j+2)^(q1+1)+(n-j)^(q1+1)-2*(n-j+1)^(q1+1))*(-35*y1(j)+17.5*y2(j)-0.5*y4(j));
        M6=M6+((n-j+2)^(q2+1)+(n-j)^(q2+1)-2*(n-j+1)^(q2+1))*(-64.062*y1(j)+31.9723*y2(j)+0.9706*(x4(j)-y4(j))+x1(j)*(x3(j)-39.031)+19.9723*x2(j));
        M7=M7+((n-j+2)^(q3+1)+(n-j)^(q3+1)-2*(n-j+1)^(q3+1))*(-y3(j)-2*x3(j)+x1(j)*x2(j));
        M8=M8+((n-j+2)^(q4+1)+(n-j)^(q4+1)-2*(n-j+1)^(q4+1))*(1.8303*y1(j)-0.9706*y2(j)-1.9723*y4(j)+0.9152*x1(j)+x2(j)*(x3(j)-0.9706)+2.4723*x4(j));
        
        M9=M9+((n-j+2)^(q1+1)+(n-j)^(q1+1)-2*(n-j+1)^(q1+1))*(-35*z1(j)+17.5*z2(j)-0.5*z4(j));
        M10=M10+((n-j+2)^(q2+1)+(n-j)^(q2+1)-2*(n-j+1)^(q2+1))*(-64.062*z1(j)+31.9723*z2(j)+0.9706*(x4(j)-z4(j))+x1(j)*(x3(j)-39.031)+19.9723*x2(j));
        M11=M11+((n-j+2)^(q3+1)+(n-j)^(q3+1)-2*(n-j+1)^(q3+1))*(-z3(j)-2*x3(j)+x1(j)*x2(j));
        M12=M12+((n-j+2)^(q4+1)+(n-j)^(q4+1)-2*(n-j+1)^(q4+1))*(1.8303*z1(j)-0.9706*z2(j)-1.9723*z4(j)+0.9152*x1(j)+x2(j)*(x3(j)-0.9706)+2.4723*x4(j));
        
        
        N1=N1+((n-j+1)^q1-(n-j)^q1)*(a*(x2(j)-x1(j))+x4(j));
        N2=N2+((n-j+1)^q2-(n-j)^q2)*((d-x3(j))*x1(j)+c*x2(j));
        N3=N3+((n-j+1)^q3-(n-j)^q3)*(x1(j)*x2(j)-x3(j)*b);  
        N4=N4+((n-j+1)^q4-(n-j)^q4)*(x3(j)*x2(j)+x4(j)*r); 
        
        N5=N5+((n-j+1)^q1-(n-j)^q1)*(-35*y1(j)+17.5*y2(j)-0.5*y4(j));
        N6=N6+((n-j+1)^q2-(n-j)^q2)*(-64.062*y1(j)+31.9723*y2(j)+0.9706*(x4(j)-y4(j))+x1(j)*(x3(j)-39.031)+19.9723*x2(j));
        N7=N7+((n-j+1)^q3-(n-j)^q3)*(-y3(j)-2*x3(j)+x1(j)*x2(j));  
        N8=N8+((n-j+1)^q4-(n-j)^q4)*(1.8303*y1(j)-0.9706*y2(j)-1.9723*y4(j)+0.9152*x1(j)+x2(j)*(x3(j)-0.9706)+2.4723*x4(j));
        
        N9=N9+((n-j+1)^q1-(n-j)^q1)*(-35*z1(j)+17.5*z2(j)-0.5*z4(j));
        N10=N10+((n-j+1)^q2-(n-j)^q2)*(-64.062*z1(j)+31.9723*z2(j)+0.9706*(x4(j)-z4(j))+x1(j)*(x3(j)-39.031)+19.9723*x2(j));
        N11=N11+((n-j+1)^q3-(n-j)^q3)*(-z3(j)-2*x3(j)+x1(j)*x2(j));  
        N12=N12+((n-j+1)^q4-(n-j)^q4)*(1.8303*z1(j)-0.9706*z2(j)-1.9723*z4(j)+0.9152*x1(j)+x2(j)*(x3(j)-0.9706)+2.4723*x4(j));
    end
    
    xx1(n+1)=x10+h^q1*N1/(q1*gamma(q1));
    xx2(n+1)=x20+h^q2*N2/(q2*gamma(q2));
    xx3(n+1)=x30+h^q3*N3/(q3*gamma(q3));
    xx4(n+1)=x40+h^q4*N4/(q4*gamma(q4));
    
    yy1(n+1)=y10+h^q1*N5/(q1*gamma(q1));
    yy2(n+1)=y20+h^q2*N6/(q2*gamma(q2));
    yy3(n+1)=y30+h^q3*N7/(q3*gamma(q3));
    yy4(n+1)=y40+h^q4*N8/(q4*gamma(q4));
    
    zz1(n+1)=z10+h^q1*N9/(q1*gamma(q1));
    zz2(n+1)=z20+h^q2*N10/(q2*gamma(q2));
    zz3(n+1)=z30+h^q3*N11/(q3*gamma(q3));
    zz4(n+1)=z40+h^q4*N12/(q4*gamma(q4));
    
    x1(n+1)=x10+h^q1*(a*(xx2(n+1)-xx1(n+1))+xx4(n+1)+M1)/gamma(q1+2);
    x2(n+1)=x20+h^q2*((d-xx3(n+1))*xx1(n+1)+c*xx2(n+1)+M2)/gamma(q2+2);
    x3(n+1)=x30+h^q3*(xx1(n+1)*xx2(n+1)-xx3(n+1)*b+M3)/gamma(q3+2); 
    x4(n+1)=x40+h^q4*(xx3(n+1)*xx2(n+1)+xx4(n+1)*r+M4)/gamma(q4+2); 
    
    y1(n+1)=y10+h^q1*(-35*yy1(n+1)+17.5*yy2(n+1)-0.5*yy4(n+1)+M5)/gamma(q1+2);
    y2(n+1)=y20+h^q2*(-64.062*yy1(n+1)+31.9723*yy2(n+1)+0.9706*(xx4(n+1)-yy4(n+1))+xx1(n+1)*(xx3(n+1)-39.031)+19.9723*xx2(n+1)+M6)/gamma(q2+2);
    y3(n+1)=y30+h^q3*(-yy3(n+1)-2*xx3(n+1)+xx1(n+1)*xx2(n+1)+M7)/gamma(q3+2); 
    y4(n+1)=y40+h^q4*(1.8303*yy1(n+1)-0.9706*yy2(n+1)-1.9723*yy4(n+1)+0.9152*xx1(n+1)+xx2(n+1)*(xx3(n+1)-0.9706)+2.4723*xx4(n+1)+M8)/gamma(q4+2); 
    
    z1(n+1)=y10+h^q1*(-35*zz1(n+1)+17.5*zz2(n+1)-0.5*zz4(n+1)+M9)/gamma(q1+2);
    z2(n+1)=y20+h^q2*(-64.062*zz1(n+1)+31.9723*zz2(n+1)+0.9706*(xx4(n+1)-zz4(n+1))+xx1(n+1)*(xx3(n+1)-39.031)+19.9723*xx2(n+1)+M10)/gamma(q2+2);
    z3(n+1)=y30+h^q3*(-zz3(n+1)-2*xx3(n+1)+xx1(n+1)*xx2(n+1)+M11)/gamma(q3+2); 
    z4(n+1)=y40+h^q4*(1.8303*zz1(n+1)-0.9706*zz2(n+1)-1.9723*zz4(n+1)+0.9152*xx1(n+1)+xx2(n+1)*(xx3(n+1)-0.9706)+2.4723*xx4(n+1)+M12)/gamma(q4+2); 
       
end
x1=x1';x2=x2';x3=x3';x4=x4';
X=[x1,x2,x3,x4];
y1=y1';y2=y2';y3=y3';y4=y4';
Y=[y1,y2,y3,y4];
z1=z1';z2=z2';z3=z3';z4=z4';
Z=[z1,z2,z3,z4];
e1=z1-y1;e2=z2-y2;e3=z3-y3;e4=z4-y4;

for n=1:(N+1)
    r(n)=sqrt((z1(n)-y1(n))^2+(z2(n)-y2(n))^2+(z3(n)-y3(n))^2+(z4(n)-y4(n))^2);
end


plot3(x3(100:end),x1(100:end),x2(100:end),'k');
xlabel('\itz','fontsize',18);
ylabel('\itx','fontsize',18);
zlabel('\ity','fontsize',18);
figure
plot(x1(100:end),x4(100:end),'k');
xlabel('\itx','fontsize',18);
ylabel('\itw','fontsize',18);
figure
plot(X(1000:end,1),Y(1000:end,1),'k')
xlabel('\itx_{1}','fontsize',18);
ylabel('\ity_{1}','fontsize',18);
figure
plot(X(1000:end,2),Y(1000:end,2),'k')
xlabel('\itx_{2}','fontsize',18);
ylabel('\ity_{2}','fontsize',18);
figure
plot(X(1000:end,3),Y(1000:end,3),'k')   
xlabel('\itx_{3}','fontsize',18);
ylabel('\ity_{3}','fontsize',18);
figure
plot(X(1000:end,4),Y(1000:end,4),'k')   
xlabel('\itx_{4}','fontsize',18);
ylabel('\ity_{4}','fontsize',18);
figure
plot(t,r,'k')
xlabel('\itt /s','fontsize',18);
ylabel('\itr','fontsize',18);
figure
plot(t,e1,'k',t,e2,'k--',t,e3,'k:',t,e4,'k-.')
xlabel('\itt /s','fontsize',18);
ylabel('\ite_{1},e_{2},e_{3},e_{4}','fontsize',18);
legend('\ite_{1}','\ite_{2}','\ite_{3}','\ite_{4}','fontsize',18)