function orbit_diagram(rmin,rmax,Ntosettle,Nrecord)
% ORBIT_DIAGRAM: Create an orbit diagram for the logistic map.

% ORBIT_DIAGRAM(rmin,rmax,Ntosettle,Nrecord) creates a orbit map analagous
% to figure 10.2.7 in Strogatz. You specify the minimum and maximum values
% of 'r' that you'd like to see, and it picks 301 intermediate values of r
% (inclusive).  It randomly selects a point in [0,1] as the initial
% condition for a given r, allows Ntosettle iterations to eliminate
% transient behavior, and then records the next Nrecord iterations for
% output.  The results are then plotted.

clf; hold on;% Clear any other figures that might be around and prepare to plot multiple points on the same figure.
r = linspace(rmin,rmax,301); % Selects all of the values of r to be examined.

for j = 1:301;
    xold = rand(1);
    for i = 1:Ntosettle;
        xnew = mapping(r(j),xold); % This loop eliminates transients by running the map Ntosettle times.
        xold = xnew;
    end

    x = zeros(1,Nrecord+1);% Initialize our results vector x by setting it to a bunch of zeros.  
    % This isn't necessary, but it allows Matlab to run quicker.
    x(1) = xold;

    for i = 2:Nrecord+1;
        x(i) = mapping(r(j),x(i-1));% Record the result of 'Nrecord' iterations.
    end

    rforplotting = r(j)*ones(1,Nrecord+1); % The abscissa vector for this value of r.
    plot(rforplotting,x,'.b','MarkerSize',5) 
end
axis([rmin rmax 0 1])
xlabel('r','FontSize',14,'FontName','Times')% A series of labelling commands, with some options changed to improve legibility.
ylabel('x','FontSize',14,'FontName','Times')
title('Orbit map for the logistic function','FontSize',16,'FontName','Times')
set(gca,'FontSize',12,'FontName','Times')% gca is the current figure's axis settings.  This changes the labels on the tick marks.

%---------------------------------------------------------------
function out = mapping(R,z);% A subfunction to define the logistic map.
out = R.*z.*(1-z);