function xlist = draw_cobweb(x0,N,xp)

% function xlist = draw_cobweb(x0,N,xp)
%
% NONLINEAR PHENOMENA: ECE 183 / ME 169 / PHYS 106
% For Recitation 1. Poincare maps and cobwebs.
% - Katie Byl, UCSB, March 28, 2013
%
% Note: Define the function "fn" within this code to match 
% any return map function you wish.  Three examples from 
% Recitation 1 are given.
%
% Inputs:
%   x0:  A real-valued initial value for the cobweb plot.
%   N :  Number of steps to simulate and draw the cobweb. (Optional.)
%   xp:  A vector of real values, for which the return map will be drawn.
%        (Optional.)
%
% Output:
%   xlist:  The sequential list of "x(n)" values, given xlist(1)=x0.

if ~exist('N','var')
    N=20;  % Number of steps to simulate.
end
if ~exist('xp','var')
    xp=0:.01:10;  % A range over which to plot RETURN MAP and set axes.
end

xp2 = fn(xp);

figure(1); clf
xmin = min(xp); xmax = max(xp);
plot([xmin xmax],[xmin xmax],'r--','LineWidth',1); hold on
axis equal
axis([xmin xmax xmin xmax])
plot(xp,xp2,'b-','LineWidth',1);

x1 = x0;
x2 = fn(x1);
plot(x1+[0 0],[0 x2],'k');
plot([x1 x2],[x2 x2],'k-');
xlist = zeros(1,N+1);
xlist(1)=x1;
xlist(2)=x2;
for n = 1:N-1
    x1 = x2;
    x2 = fn(x1);
    xlist(2+n)=x2;
    plot([x1 x1],[x1 x2],'k');
    plot([x1 x2],[x2 x2],'k');
    
end
end

function x2 = fn(x1)

% Recitation 1: First (piece-wise) return map.
x2 = (x1<4).*(-1+1.5*x1) + (x1>=4).*(5-1.5*(x1-4));

% Recitation 1: Second return maps.
% x2 = 2*atan(x1);
% x2 = -2*atan(x1);

end





