function lorenzsys
  % Three-dimensional plot of solution of Lorenz equations
  rho   = 28;
  sigma = 10;
  beta  = 8/3;
  tspan = linspace(0, 100, 1e4 );
  y0    = [ 27, -8, 8 ];

  [ tplot, yplot ] = ode45( @(t,y)lorenzeqs(t,y,rho,sigma,beta), tspan, y0 );
  % first few points of the time history
  plot( tplot(1:200), yplot(1:200,1),'k-', tplot(1:200), yplot(1:200,2), 'k--', tplot(1:200), yplot(1:200,3), 'k-.' )
  set(gca,'fontsize',16);
  figure
  % phase diagram
  plot3( yplot(:,1), yplot(:,2), yplot(:,3), '-k' )
  set(gca,'fontsize',16);


function ydot=lorenzeqs(t,y,rho,sigma,beta)
  ydot(1,:) = -beta*y(1,:) + y(2,:).*y(3,:);
  ydot(2,:) = sigma*( y(3,:)-y(2,:) );
  ydot(3,:) = -y(2,:).*y(1,:) + rho*y(2,:)-y(3,:);
end

end