function cobweb(fun,x0,xmin,xmax,n,tol)
% cobweb(fun,x0,xmin,xmax,n,tol)
%
% COBWEB plots the orbit of X0 for any given iteration function FUN.  
% XMIN and XMAX determine the range of the plot; default range is [0 1.5]. 
% N is the maximum number of iterations to do; default 30.  TOL  
% determines how close to a fixed point you must get to exit (XMIN and 
% XMAX must be given in order to also specify N and TOL); default TOL is
% 1.e-3.   
% 
% Typing COBWEB by itself performs the following example:  
%     cobweb('2.8*x*(1-x)',.5,0,1.5,30,.001)
% Other examples:
%     cobweb('5+x-x^2',2.5,0,5,3)
%     cobweb('5/x',2.5,1.5,3,5)
%     cobweb('1+x-x^2/5',2.5,2.2,2.5,10,1e-6)
%     cobweb('(x+5/x)/2',2.5,2.2,2.5,10,1e-6)
  
  if nargin < 1
    fun = '2.8*x*(1-x)';
  end
  
  if nargin < 2
    x0= .5;
  end
  
  if nargin < 3
    xmin = 0;
    xmax = 1.5;  
  end
  
  if nargin < 5
    n=30;
  end
  
  if nargin < 6
    tol = 1.e-3;
  end
  
  clf
  ezplot('x', [xmin xmax]);
  hold on;
  axis equal
  ezplot(fun, [xmin xmax]);
  axis([xmin xmax xmin xmax]);
  
  x = x0;
  plot(x,max(0,xmin),'ro')
  pause
  plot(x,max(0,xmin),'go')
  
  [F,msg] = fcnchk(fun,'vectorized'); % Make a vectorized feval-able function.
  if ~isempty(msg), error(msg); end
  
  mes='';
  % initially start vert line at y = 0:
  y = F(x);
  plot([x x],[0 y],'--r')
  plot([min(x,y) max(x,y)],[y y],'--r')
  plot(y,max(0,xmin),'ro')
  pause
  plot([x x],[0 y],'--g')
  plot([min(x,y) max(x,y)],[y y],'--g')
  plot(y,max(0,xmin),'go')
  
  for i = 1:n-1,
    xold = x;
    x = y;
    y = F(x);
    plot([x x],[min(x,y) max(x,y)],'--r')  % vert line
    plot([min(x,y) max(x,y)],[y y],'--r')  % horz line  
    plot(y,max(0,xmin),'ro')
    pause
    plot([x x],[min(x,y) max(x,y)],'--g')  % vert line
    plot([min(x,y) max(x,y)],[y y],'--g')  % horz line  
    plot(y,max(0,xmin),'go')
    if abs(x-y) <= tol
      mes=sprintf('\nFixed point reached at x=%g at %d iteration(s).',x,i+1);
      disp(mes);
      break 
    elseif abs(x-y)>abs(x-xold)
      mes=sprintf('\nIterations diverged.');
      disp(mes);
      break 
    end
  end
  
  
  if i == n-1
    mes=sprintf('\nMaximum iterations exceeded.');
    disp(mes);
  end
  
  str = sprintf('Orbit of x=%g in f(x) = %s%s',x0,fun,mes);
  title(str);
  xlabel('x');
  ylabel('y');
  hold off;