clear all;
close all;
x(1,:)=[3 0.5 8 7];
h=0.05;
for n=1:10000
    if sqrt(x(n,1)^2+x(n,2)^2+x(n,3)^2+x(n,4)^2)<=8
        a=-6;%不穩定的線性系統，系統Jacobian矩陣的特征值為：-3，-2，1，1）
        b=7;
        c=3;
        d=-3;
        
    K1=x(n,2);
    L1=x(n,3);
    M1=x(n,4);
    N1=a*x(n,1)+b*x(n,2)+c*x(n,3)+d*x(n,4);
    
    K2=x(n,2)+h*L1/2;
    L2=x(n,3)+h*M1/2;
    M2=x(n,4)+h*N1/2;
    N2=a*(x(n,1)+h*K1/2)+b*(x(n,2)+h*L1/2)+c*(x(n,3)+h*M1/2)+d*(x(n,4)+h*N1/2);
    
    K3=x(n,2)+h*L2/2;
    L3=x(n,3)+h*M2/2;
    M3=x(n,4)+h*N2/2;
    N3=a*(x(n,1)+h*K2/2)+b*(x(n,2)+h*L2/2)+c*(x(n,3)+h*M2/2)+d*(x(n,4)+h*N2/2);
    
    K4=x(n,2)+h*L3;
    L4=x(n,3)+h*M3;
    M4=x(n,4)+h*N3;
    N4=a*(x(n,1)+h*K3)+b*(x(n,2)+h*L3)+c*(x(n,3)+h*M3)+d*(x(n,4)+h*N3);
    
          x(n+1,:)=x(n,:)+[h*(K1+2*K2+2*K3+K4)/6 h*(L1+2*L2+2*L3+L4)/6 h*(M1+2*M2+2*M3+M4)/6 h*(N1+2*N2+2*N3+N4)/6];
    else
        a=-6;    %穩定的線性系統，系統Jacobian矩陣的特征值為：-3，-2，-1，-1）
        b=-17;
        c=-17;
        d=-3;
        
    K1=x(n,2);
    L1=x(n,3);
    M1=x(n,4);
    N1=a*x(n,1)+b*x(n,2)+c*x(n,3)+d*x(n,4);
    
    K2=x(n,2)+h*L1/2;
    L2=x(n,3)+h*M1/2;
    M2=x(n,4)+h*N1/2;
    N2=a*(x(n,1)+h*K1/2)+b*(x(n,2)+h*L1/2)+c*(x(n,3)+h*M1/2)+d*(x(n,4)+h*N1/2);
    
    K3=x(n,2)+h*L2/2;
    L3=x(n,3)+h*M2/2;
    M3=x(n,4)+h*N2/2;
    N3=a*(x(n,1)+h*K2/2)+b*(x(n,2)+h*L2/2)+c*(x(n,3)+h*M2/2)+d*(x(n,4)+h*N2/2);
    
    K4=x(n,2)+h*L3;
    L4=x(n,3)+h*M3;
    M4=x(n,4)+h*N3;
    N4=a*(x(n,1)+h*K3)+b*(x(n,2)+h*L3)+c*(x(n,3)+h*M3)+d*(x(n,4)+h*N3);
    
      x(n+1,:)=x(n,:)+[h*(K1+2*K2+2*K3+K4)/6 h*(L1+2*L2+2*L3+L4)/6 h*(M1+2*M2+2*M3+M4)/6 h*(N1+2*N2+2*N3+N4)/6];
  end
  if n>=1000
      y(n,:)=x(n,:);
  end
  n
end
plot3(y(:,1),y(:,2),y(:,3),'k') ;%系統的四個Lyapunov指數為：1.8074,0.0507,-0.8796,-3.5796.  
xlabel('x1');
ylabel('x2');
zlabel('x3');