% Numerical solutions for the Lorenz equations

% This is a script that we went over in class on 3/26. To run all sections 
% at once, click 'Run'. Otherwise, click 'Run and Advance'.

%% Run 1

y0 = [.2 .3 .4]; 
options = odeset('RelTol',.0000001,'AbsTol',.000000001);
r = 1.5;
tspan = 0:.01:50;
[tt,yy] = ode45('lorenz_solver',tspan,y0,options,r);
figure, plot3(yy(:,1),yy(:,2),yy(:,3));

%% Run 2

y0 = [1 1 1]; 
options = odeset('RelTol',.0000001,'AbsTol',.000000001);
r = 8;
tspan = 0:.01:50;
[tt,yy] = ode45('lorenz_solver',tspan,y0,options,r);
figure, plot3(yy(:,1),yy(:,2),yy(:,3));

%% Run 3

y0 = [6.7 6.7 17]; 
options = odeset('RelTol',.0000001,'AbsTol',.000000001);
r = 18;
tspan = 0:.01:50;
[tt,yy] = ode45('lorenz_solver',tspan,y0,options,r);
figure, plot3(yy(:,1),yy(:,2),yy(:,3));

%% Run 4

y0 = [9 8 27]; 
options = odeset('RelTol',.0000001,'AbsTol',.000000001);
r = 28;
tspan = 0:.01:200;
[tt,yy] = ode45('lorenz_solver',tspan,y0,options,r);
figure, plot3(yy(:,1),yy(:,2),yy(:,3));

%% Run 5

% Initial condition
y0 = [9 8 27]; 
options = odeset('RelTol',.0000001,'AbsTol',.000000001);
r1 = 28;
tspan = 0:.01:50;
[tt1,yy] = ode45('lorenz_solver',tspan,y0,options,r1);

% Small perturbation of our initial condition 
z0 = [9 8 27.0001]; 
options = odeset('RelTol',.0000001,'AbsTol',.000000001);
r2 = 28;
tspan = 0:.01:50;
[tt2,zz] = ode45('lorenz_solver',tspan,z0,options,r2);
figure, hold on
plot3(yy(:,1),yy(:,2),yy(:,3),'r');
plot3(zz(:,1),zz(:,2),zz(:,3),'b');

%% Run 6
% Try several values of tmax, to see what happens/when the solutions
% diverge. You could also try varying the perturbation size from the
% initial condition. 'Sensitive dependence on initial conditions,' one of 
% the hallmarks of chaotic behavior in continuous systems, depends on both.
% Some tmax values to try for the different sets of initial conditions:
% tmax = 5, 10, 20, 30, ... 35, 40

tmax = 5;
offset = 0.01;

y0 = [9 8 27]; 
options = odeset('RelTol',.0000001,'AbsTol',.000000001);
r1 = 28;
tspan = 0:.01:tmax;
[tt1,yy] = ode45('lorenz_solver',tspan,y0,options,r1);

znew = 27+offset;
z0 = [9 8 znew];
options = odeset('RelTol',.0000001,'AbsTol',.000000001);
r2 = 28;
tspan = 0:.01:tmax;
[tt2,zz] = ode45('lorenz_solver',tspan,z0,options,r2);
figure, hold on
plot3(yy(:,1),yy(:,2),yy(:,3),'r');
plot3(zz(:,1),zz(:,2),zz(:,3),'b');