% Logistic map

%% Taken from http://www.physics.nau.edu/~hart/classes/600_fall_2004/matlab/logmap_orbit_plot.m

clc; clear; close

N = 400;  % Number of iterations for each value of r
M = 300;  % Number of different values of r
np = 200; % Number of final iterates to plot for each value of r 

r_i = 3.4;  % Starting value of r
r_f =   4;  % Ending value of r
yll = -.1;  % Lower y limit of plot
yul =   1;  % Upper y limit of plot

%%% Zoom
r_i = 3.847;  % Starting value of r
r_f = 3.857;  % Ending value of r
yll = .13;    % Lower y limit of plot
yul = .185;   % Upper y limit of plot

plotting_x = r_i;
plotting_y = 0;
fh=plot([r_i,0],'.','EraseMode','none'); % Make the initial plot and get the figure handle
axis([r_i,r_f, yll, yul])                % Set the axis limits for the plot

r=linspace(r_i,r_f,M);                % Create an array of M values of r equally spaced r_i < r(n) < r_j

for j = 1:M                           % Loop over each value of r
    x = zeros(1,N);                   % Create an array of x values (much faster than appending as we go
    x(1) = .6;                        % Starting value of x(n) is arbitrary
    for n = 1:N-1                     % For each value of r, iterate x long enough for the transient to die away
        x(n+1) = r(j)*x(n)*(1-x(n));  % Get next value of x(n+1). See Eq. 1 in Chap. 10 of Strogatz
    end
    rdata = ones(1,np).*r(j);
    fx = x((N-np)+1:N);
    plotting_x = [plotting_x, rdata];
    plotting_y = [plotting_y, fx];
    set(fh,'xdata',rdata,'ydata',fx); % Add the new points to the plots
    drawnow                           % This command causes Matlab to immediately draw the new points
end
clf;
plot(plotting_x,plotting_y,'.')
axis([r_i,r_f, yll, yul])
