function doall=lorenz
clear all; close all;
% Eli Tziperman: example for course in nonlinear dynamics and
% chaos.
% + additional changes by Vered to get the unstable manifolds
% Solve Lorenz system using ode45 (Runga-Kutta package of matlab) 

% the equation:
%  dx/dt = sigma*(y-x)
%  dy/dt = r*x-y-x*z
%  dz/dt = x*y-b*z

% define parameters:
% ------------------
sigma=10.0;
b=8.0/3.0;
r=15;

delta_t=0.01;

% integration time:
T=5.0;
%set error bounds for integration
options=odeset('RelTol',10^(-6),'AbsTol',[10^(-6)*ones(1,3)]);


% integrate:
% initial conditions (x,y,z)=[0.00001 0 0]
 [t,x] = ode45(@lor,[0 T],[0.00001 0 0],options,[sigma; b ;r]);
 % initial conditions (x,y,z)=[-0.00001 0 0]
 [t,xl] = ode45(@lor,[0 T],[-0.00001 0 0],options,[sigma; b ;r]);

% plot simple time series:
plot(t(:),x(:,1))
plot(t(:),xl(:,1),'r')
xlabel('time')
ylabel('x')

% plot phase space:
figure
plot(x(:,1),x(:,3))
hold
plot(xl(:,1),xl(:,3),'r')
xlabel('x')
ylabel('z')

% plot phase space:
figure
plot(x(:,1),x(:,2))
hold
plot(xl(:,1),xl(:,2),'r')
xlabel('x')
ylabel('y')

% plot phase space:
figure
plot3(x(:,1),x(:,2),x(:,3))
hold
plot3(xl(:,1),xl(:,2),xl(:,3),'r')
xlabel('x')
ylabel('y')
zlabel('z')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy = lor(t,x,kh)
sigma=kh(1); b=kh(2); r=kh(3);
dy=[sigma*(x(2)-x(1)),r*x(1)-x(2)-x(1)*x(3),x(1)*x(2)-b*x(3)]';