% logistic map: x(n+1) = A * x(n) * ( 1 - x(n))
% ========================================================================
% Nonlinear Dynamical Systems
% date : January 23,2006, Monday
% IES: http://www.cavs.msstate.edu/hse/ies/projects/nsf_nonlinear
% Created by: Sundararajan Srinivasan and Saurabh Prasad
% ========================================================================
clear;
clc;

a = 4;
x_axis = [0 : 0.01: 1.0];
y_axis = x_axis;
plot(x_axis, y_axis);
hold on;
plot(x_axis, a*x_axis.*(1-x_axis));
hold on;

N = 100;
N0 = 100;
num_steps = 40;

x0 = 0.1;
x= x0;
y = 0;
for i = 1 : N
    y_next = a * x * (1 - x);
    plot([x x], [y y_next]);
    hold on;
    x_next = y_next;
    plot([x x_next], [y_next y_next]);
    hold on;
    
    x = x_next;
    y = y_next;
end; 
    
figure;
subplot(2,1,1);
hold on;
lyp =[];

a = 2.8 : 0.005 : 4;
for i = 1:length(a)
    x0 = 0.1;
    x= x0;
    y = 0;
    x_accum = [];
    x_accum1 = [];
    
    for j = 1 : N + N0
        y_next = a(i) * x * (1 - x);
        x_next = y_next;
        x = x_next;
        y = y_next;
        
        x_accum1 = [x_accum1 x];
        if j > N0
            x_accum = [x_accum x];
        end;
    end;
    A = ones(size(x_accum)) * a(i);
    plot(A, x_accum, '.');
    
    lyp = [lyp lyapunov(x_accum1(:), 3, 2, 1, num_steps)];
end;
subplot(2,1,2);
hold on;
plot(a, zeros(size(a)));
plot(a, lyp)
