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

clear all
close all

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

%parameters
a = 8;
b = 1;
p = [a;b];

%Equilibria
xeq = a/5;
yeq = 1+xeq^2;
Ueq = [xeq;yeq];

% [X,Y] = meshgrid([0:0.1:10],[0:0.1:10]);
% F = bz(0,[X,Y],p);
% 
% figure(4)
% quiver(X,Y,F1,F2)



%Initial condition(s)
U1 =[0.5,0];
U2 =[0,3]; 
%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) bz(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(xeq,yeq,'.k','MarkerSize',10)
hold off
xlabel('x')
ylabel('y')