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



%Run Time
tf = 100;
%tf=17;
%tf =15;

%parameters
c = 2.6;
p = c;


x = linspace(-1,2,100);
y = linspace(-1,1,100);
[X,Y] = meshgrid(x,y);
Z = Y.^2/2 + X.^2/2 - X.^3/3;
figure(3)
contour(X,Y,Z,40)


lam_u = (-c +sqrt(c^2 +4))/2;
V_u = [1; lam_u];
X1 = [1;0];

%Initial condition(s)
U1 = X1 + 0.01*V_u;%[-0.4165 -0.9091] %0];%[0 1 0]
U2 = X1-0.01*V_u;%[0 -1 0];  %%Trajectories on both sides of the unstable manifold Ws(C0)


% Generate Solution


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

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

% Create 2d Plot of solutions
figure(1)
plot(y1(:,1),y1(:,2),'linewidth',3)
hold on
plot(y2(:,1),y2(:,2),'linewidth',3)
plot([1 0],[0 0],'.k','MarkerSize',15)
hold off
xlim([-0.4 1.5])
ylim([-0.5 0.5])