function this=void
clear all; set(0,'defaultaxesfontsize',15);format long
%%% behaviour of the chaotic (deterministic) Lorenz ' 63 system
%%% based on Ex 2.6 in: Law, Stuart, Zygalakis (2015), "Data Assimilation"

%set up the parameters of the problem
a=10;
b=8/3;
r=28;
sd=1; rng(sd); %always using the same random numbers
initial=randn(3,1);
initial1=initial;
initial1(1)=initial(1)+0.0001; %perturbed initial condition

%calculate the trajectory 
[t1,y]=ode45(@(t,y) lorenz63(t,y,a,b,r), [0 100], initial);
[t,y1]=ode45(@(t,y) lorenz63(t,y,a,b,r), t1, initial1);

for i=1:length(t)
error(i)=norm(y(i,:)-y1(i,:),2);
end

figure(1) 
semilogy(t,error,'b'); axis([0 100 10^-6 10^2])
title('l2-deviation between trajectories')
set(gca,'YTick',[10^-6 10^-4 10^-2 10^0 10^2])
set(gca,'XTickLabel',{'-pi','-pi/2','0','pi/2','pi'})

figure(2) 
plot3(y(:,1),y(:,2),y(:,3))
hold on
plot3(y1(:,1),y1(:,2),y1(:,3),'r')
hold off
title('strange attractor of the Lorenz 63 system')

figure(3)
plot(t,y(:,1)); xlabel('time'); title(['first component'])



function dy=lorenz63(t,y,a,b,r)

dy(1)=a*(y(2)-y(1));
dy(2)=-a*y(1)-y(2)-y(1)*y(3);
dy(3)=y(1)*y(2)-b*y(3)-b*(r+a);

dy=dy';