% ###  EXlogisticBIF.m ###        10.16.14
% Produces a bifurcation diagram for the logistic map
%   x_{n+1} = x_n * r * (1- x_n)
% by varying r over the range leading to the period-doubling cascade
% --> also visualizes the associated 'time course'

clear;  figure(1); clf; hold on;
% ---------------
% User inputs
range= [2 4];   % min and max values to compute bifurcation diagram over [2 4]
Nr= 200;     % # of steps over range [100]
x0= 0.1;    % starting x value [0.1]
Nsettl= 50; % # of runs allowed for 'settling' [50]
Nit= 100;   % # of iterations to plot for a given value of r [200]
rPlot= 3.9; % for 'timecourse' plot, specify associated r value (must be inside range!)
% ---------------
rmin= range(1); rmax= range(2);
% loop through each r value
for nn=1:Nr
    r(nn)= rmin + nn*(rmax-rmin)/Nr;    % update r
    x= x0;  % reset to IC
    xS(1)= x;   % store first point
    indx=2; % reset indexer (for 2nd iterate) 
    for mm=0:Nsettl+Nit     % loop through the iterations of the map
        x= r(nn)*x*(1-x);       % deal with mapping
        xS(nn,indx)= x;        % store values
        indx= indx+1;       % update indexer
    end
    % plot points for a given iteration *past* the settling time
    plot(r(nn)*ones(Nit+1),xS(nn,Nsettl:Nsettl+Nit),'k.')  
end
% ----
xlabel('r');    ylabel('x_n')
title('Bifurcation Diagram for the Lositic Map [x_{n+1} = r*x_n*(1-x_n)]')
% ----
% also plot x_n as function of n for relevant r value (as specified)
[junk indxR] = min(abs(r-rPlot));   % search for closest r value to rPlot
n= linspace(0,size(xS,2),size(xS,2));
figure(2); clf;
plot(n,xS(indxR,:),'kd-'); hold on;
xlabel('n');    ylabel('x_n');
stem(Nsettl,max(xS(indxR,:)),'r-','marker','none'); % indicate bound for 'settling'
figure(1); stem(rPlot,max(xS(:)),'r-','marker','none'); % indicate r for which 'time course' is plotted