a=0.2;
b=0.2;
c=4.17;     % 2.5 (1) 3.5 (2) 4 (4) 4.15 (8) 4.17 (16) 5 (chaos)
dt=0.01;
t0=600; if(c>4.3) t0=t0; end

clf
drawnow;

kx=zeros(4,1);
ky=zeros(4,1);
kz=zeros(4,1);
psx=zeros(0);
psz=zeros(0);

% Start at the origin
x=0; x1=x;
y=0; y1=y;
z=0; z1=z;

X=[];
Y=[];
Z=[];
xm=[];


% Integrate forward in time
n=0;
t=0;
while(t<t0)
    n=n+1;
    t=t+dt;

    x2=x1;
    y2=y1;
    z2=z1;

    x1=x;
    y1=y;
    z1=z;

% do a Runge-Kutta integration step
    v(1)=-y-z;
    v(2)=x+a*y;
    v(3)=b+z*(x-c);
    kx(1)=dt*v(1);
    ky(1)=dt*v(2);
    kz(1)=dt*v(3);
    x=x1+0.5*kx(1);
    y=y1+0.5*ky(1);
    z=z1+0.5*kz(1);
%
    v(1)=-y-z;
    v(2)=x+a*y;
    v(3)=b+z*(x-c);
    kx(2)=dt*v(1);
    ky(2)=dt*v(2);
    kz(2)=dt*v(3);
    x=x1+0.5*kx(2);
    y=y1+0.5*ky(2);
    z=z1+0.5*kz(2);
%
    v(1)=-y-z;
    v(2)=x+a*y;
    v(3)=b+z*(x-c);
    kx(3)=dt*v(1);
    ky(3)=dt*v(2);
    kz(3)=dt*v(3);
    x=x1+kx(3);
    y=y1+ky(3);
    z=z1+kz(3);
%
    v(1)=-y-z;
    v(2)=x+a*y;
    v(3)=b+z*(x-c);
    kx(4)=dt*v(1);
    ky(4)=dt*v(2);
    kz(4)=dt*v(3);
    x=x1+(kx(1)+2*kx(2)+2*kx(3)+kx(4))/6;
    y=y1+(ky(1)+2*ky(2)+2*ky(3)+ky(4))/6;
    z=z1+(kz(1)+2*kz(2)+2*kz(3)+kz(4))/6;

% Find successive maxima of x(t)
    if(x1>x & x1>x2)
      xm=[xm x];
      subplot(2,2,4);
      hold on;
      plot([0 12],[0 12],'r');
      m=length(xm);
      if(m>1)
        plot(xm(m-1),xm(m),'b.','markersize',1);
      end
      title('Map of consecutive maxima');
      xlabel('x_n');
      ylabel('x_n_+_1');
      drawnow
    end

% Find intersections with the Poincare surface of section (x=0)
  if(t>200)
    if(x1*x<0 & y<y1)
      psz=[psz; y z];
    end
    if(x1*x<0 & y>y1)
      psx=[psx; y z];
      subplot(2,2,3);
      hold off;
      mx=length(psx(:,1));
      mz=length(psx(:,1));
      plot([2 8],[-2 8],'w.','markersize',1);
      hold on;
      if(mx>1 & mz>1)
        plot(psx(:,1),psx(:,2),'b.','markersize',1);
        plot(psz(:,1),psz(:,2),'r.','markersize',1);
      end
      if(mx>1)
        plot(psx(mx-1,1),psx(mx-1,2),'gx');
        plot(psx(mx,1),psx(mx,2),'rx');
      end
      title('x = 0 Poincare section');
      xlabel('y');
      ylabel('z');
      drawnow
    end

    if(mod(n,4)==0)
        X=[X; x];
        Y=[Y; y];
        Z=[Z; z];
    end

% Plot trajectory in 3D and 2D projection
    if(mod(n,2000)==0);
      subplot(2,2,1);
      plot3(X,Y,Z,'r');
      title('trajectory in 3D');
      xlabel('x');
      ylabel('y');
      zlabel('z');
      drawnow
%
      subplot(2,2,2);
      hold off
      plot([0 10],[0 10],'w.','markersize',1);
      hold on
      plot(X,Y,'b');
      title('x-y projection');
      xlabel('x');
      ylabel('y');
      drawnow;
%
    end
  end
end
drawnow