a=0.2;
b=0.2;
c=5;        % 2.5 (small period limit cycle), 5 (chaos)
dt=0.01;
N=200;
l=1e-1;
t0=100; if(c>4.3) t0=1*t0; end

% Define an array of points to serve as initial conditions
X=linspace(0,l,N);
Y=0*X;
Z=0*X;

clf
drawnow;

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

XX=[];
YY=[];
ZZ=[];

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

    x=X(i);
    y=Y(i);
    z=Z(i);
    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(i)=x1+(kx(1)+2*kx(2)+2*kx(3)+kx(4))/6;
    Y(i)=y1+(ky(1)+2*ky(2)+2*ky(3)+ky(4))/6;
    Z(i)=z1+(kz(1)+2*kz(2)+2*kz(3)+kz(4))/6;
  end

  XX=[XX X(1)];
  YY=[YY Y(1)];
  ZZ=[ZZ Z(1)];
  
  if(mod(n,20)==0);
% Plot trajectory starting at the origin (blue) and stretched curve (red) 
      subplot(1,2,1);
      hold off
      plot3(X,Y,Z,'r','linewidth',2);
      hold on
      plot3(XX,YY,ZZ,'b');
      title('stretched curve');
      xlabel('x');
      ylabel('y');
      zlabel('z');
      drawnow

% Compute the length of the stretched curve
      l=0;
      for i=2:N
        l=l+sqrt((X(i)-X(i-1))^2+(Y(i)-Y(i-1))^2+(Z(i)-Z(i-1))^2);
      end
      subplot(1,2,2);
      semilogy(t,l,'r*','markersize',4);
      hold on
      title('curve length');
      xlabel('time');
      ylabel('length');
      drawnow

  end
end
drawnow