function lorenz2d
  gamma=10;
  b=8/3;
  mu=28;
  nstep=10000;
  dt=0.01;
  
  x=zeros(nstep+1);
  y=zeros(nstep+1);
  z=zeros(nstep+1);
  
  x(1)=0;
  y(1)=1.001;
  z(1)=0;

  for n=2:nstep+1
    x(n)=x(n-1)-dt*gamma*(x(n-1)-y(n-1));
    y(n)=y(n-1)+dt*(mu*x(n-1)-y(n-1)-x(n-1)*z(n-1));
    z(n)=z(n-1)+dt*(-b*z(n-1)+x(n-1)*y(n-1));
  end

  plot(x,y,x(nstep),y(nstep),'mo','markersize',10)

