function ShowLorenz()
% Framework for studying the Lorenz system
%
% [ y1'(t) ] [ beta 0 y2(t) ] [ y1(t) ]
% [ y2'(t) ] = [ 0 -sigma sigma ] [ y2(t) ]
% [ y3'(t) ] [ -y2(t) rho -1 ] [ y3(t) ]

% where the three parameters are given by (for instance)
sigma = 10;
rho = 28;
beta = 8/3;

% Note that if
eta = sqrt(beta*(rho-1));
% and
Attractor_1 = [rho-1;eta;eta];
Attractor_2 = [rho-1;-eta;-eta];
% then

% [ beta 0 y2(t) ] [ 0 ]
% [ 0 -sigma sigma ] * Attractor_i = [ 0 ]
% [ -y2(t) rho -1 ] [ 0 ]

% Integrate from 0 to tFinal with initial condition y0. Here are
% some illustrative choices

y0 = Attractor_1 + [0; 0; 3]; tFinal = 50;
% y0 = Attractor_2 + [0; 0; 3]; tFinal = 50;
% y0 = Attractor_2 + [0; 0; .3]; tFinal = 50;
% y0 = (Attractor_1 + Attractor_2)/2 + .000001*randn(3,1); tFinal = 50;

close all

% We will use ode45...
opts = odeset('RelTol',.000001);
[tVals,yVals] = ode45(@(t,y) Lorenz(t,y,beta,sigma,rho),[0 tFinal],y0,opts);
% A dynamic display...
comet3(yVals(:,1),yVals(:,2),yVals(:,3))
shg

function yp = Lorenz(t,y,beta,sigma,rho)
yp = [-beta*y(1)+y(2)*y(3) ;...
-sigma*y(2)+sigma*y(3) ; ...
-y(2)*y(1)+rho*y(2)-y(3)...
];