%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         5. practice         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Nonlinear dynamics systems  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         2019.11.13.         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       Lorenz system         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Parameters:

s=10;                       % 
r=28;                       % 
b=8/3;                      % 
start1=[5 5 5];             % 1. start point
start2=[4.9 5 5];           % 2. start point
tmax=50;                    % simulation time
tlim=[0 tmax];              % borders for plotting
xl=[-20,20];                % 
yl=[-30,30];                % 
zl=[-20,50];                % 

Lorenz_equ=@(t, xyz) [s*(xyz(2)-xyz(1));...
                           r*xyz(1)-xyz(2)-xyz(1)*xyz(3);...
                           xyz(1)*xyz(2)-b*xyz(3)]; % Lorenz equations
% Lorenz equations

figure;
    [t, x] = ode45(Lorenz_equ, tlim, start1);              % ode45 for the 1st initial position    
    [t2, x2] = ode45(Lorenz_equ, tlim, start2);           % ode45 for the 2nd initial position 

    

% plot x,y,z in the function of time (subplots)
% t-x
subplot(6,2,1);
plot(t,x(:,1), t2,x2(:,1));
xlabel('t');
ylabel('x');
xlim(tlim);
ylim(xl);

% t-y
subplot(6,2,3);
plot(t,x(:,2), t2,x2(:,2));
xlabel('t');
ylabel('y');
xlim(tlim);
ylim(yl);

% t-z
subplot(6,2,5);
plot(t,x(:,3), t2,x2(:,3));
xlabel('t');
ylabel('z');
xlim(tlim);
ylim(zl);

% x-y-z 3D plot
subplot(6,2,7:12);
plot3(x(:,1),x(:,2),x(:,3), x2(:,1),x2(:,2),x2(:,3));
xlabel('x');
ylabel('y');
zlabel('z');
xlim(xl);
ylim( yl);
zlim(zl);

subplot(6,2, [2,4,6]);
[my_what, my_where]=findpeaks(x(:,3)); %build in function from the signal processing toolbox: finds local maximums in a vector
peak_n1=my_what(1:end-1); %plot the 2 datasets in the function of each other
peak_n2=my_what(2:end);

[my_what2, my_where2]=findpeaks(x2(:,3)); %build in function from the signal processing toolbox: finds local maximums in a vector
peak2_n1=my_what2(1:end-1); %plot the 2 datasets in the function of each other
peak2_n2=my_what2(2:end);

plot(peak_n1, peak_n2, 'b.', peak2_n1, peak2_n2, 'r.'); %plot local maximums: every n+1. in the function the n.
xlim([30 50]);
ylim([30 50]);
xlabel('Zn');
ylabel('Zn+1');
title('Lorenz peak map');







