%%%%Script to solve lorenz equations,
%%%% Change constants in function file lorenz.m


%Run Time
tf = 300;
%tf=10;
%parameters
sigma = 10;
r = 28;
b = 8/3;

p = [sigma,r,b];

%%Equilibria
C0 = [0,0,0]';
Cm = [-sqrt(b*(r-1)),-sqrt(b*(r-1)),r-1]';
Cp = [+sqrt(b*(r-1)),+sqrt(b*(r-1)),r-1]';


%%Jacobian Matrix
xs = C0(1);
ys = C0(2);
J = [-sigma sigma 0 ; r -1 -xs;ys xs -b];
[VV,DD] = eig(J)

%Initial condition(s)
d0 = 6e-16;
delta = d0*[1,0,0]; %%%small pertubation of size d0
U1 = 0.2*[-0.4165 -0.9091 0]+rand(1,3)*0.00001;%[0 1 0]
U2 = U1+d0;%[0 -1 0];  %%Trajectories on one side of the unstable manifold Ws(C0)


% Generate Solution


ff =@(tt,yy) lorenz(tt,yy,p);

opts = odeset('MaxStep',2e-2)
%%Solutions for different initial conditions
[t1, y1] = ode45(ff, [0 tf], U1,opts);
[t2, y2] = ode45(ff, [0 tf], U2,opts);

% Create 3d Plot
figure(1)
plot3(y1(:,1), y1(:,2), y1(:,3));
hold on
plot3(y2(:,1), y2(:,2), y2(:,3));
plot3(C0(1),C0(2),C0(3),'.k','MarkerSize',10)
plot3(Cp(1),Cp(2),Cp(3),'.k','MarkerSize',10)
plot3(Cm(1),Cm(2),Cm(3),'.k','MarkerSize',10)
hold off
xlabel('x')
ylabel('y')
zlabel('z')




%%Plot times series in x and z
figure(2)
subplot(2,1,1)
plot(t1,y1(:,1))
xlabel('t')
ylabel('x')
subplot(2,1,2)
plot(t1,y1(:,3))
xlabel('t')
ylabel('z')


%%Plot Error
DEL = y1 - y2;
NDEL = vecnorm(DEL,2,2);
figure(3)
plot(t1,NDEL)

figure(4)
plot(t1(3000:9000),log(NDEL(3000:9000)),'LineWidth',3)
