function cobweb(a,b,x0,Niter)
% COBWEB: Create a cobweb plot for a discrete dynamical mapping function.

% COBWEB(a,b,x0,Niter) creates a cobweb plot for a mapping function 
% (usually defined in the MAPPING subfunction below).  The domain and 
% range for the plot are [a,b], and the initial condition for the system 
% is set by 'x0'.  The mapping is performed a total of 'Niter' iterations.


x=zeros(1,Niter+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)=x0; % The first entry is our given initial condition.

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

if any(x > b | x < a)
    sprintf('The maximum value of x_n is %2.4f and the minimum value is %2.4f.',[max(x) min(x)])
    error('Please adjust your domain accordingly.')
end

clf; hold on; % Clear any other figures that might be around and prepare to plot multiple curves on the same figure.

x4plotting = linspace(a,b,100); % x4plotting is a temporary series of points that allows us to...
plot(x4plotting,mapping(x4plotting),'b-')% plot the dynamical mapping function smoothly.

% Alternately, we can replace the above 2 lines with Matlab's function
% plotting command fplot:
% fplot('mapping',[a b],'b-') 

axis equal; axis([a b a b]); % Set the axes to a nice square.

plot([a b], [a b], 'k--')% Plot the line x_{n+1} = x_n for reference.

y=x(:,floor([1:.5:Niter+1])); % A handy trick to allow us to plot the cobweb. 
plot(x(1),a,'ro') % Plot the initial point as a red circle.
plot(y(1:2),[a y(3)],'r:');% Plot the initial leg with a dotted line.
plot(y(2:2*Niter),y(3:2*Niter+1),'r-.');% Plot the rest of the iterations with dot-dashed line.

xlabel('x_n','FontSize',14,'FontName','Times') % A series of labelling commands, with some options changed to improve legibility.
ylabel('x_{n+1}','FontSize',14,'FontName','Times')
title('Cobweb diagram','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(z); % A subfunction that defintes the dynamical mapping function.
out = 3.2.*z.*(1-z); % In this case, the logistic map.