function cobweb(fun,xzero,mu,nstep);

%   COBWEB(FUN,XZERO,MU,NSTEP) shows the cobweb plot of the given
%   function FUN with parameter MU, and follows a point XZERO for NSTEP 
%   mappings

%   Copyright 2003
%   Centre for Nonlinear Dynamics in Physiology and Medicine

for i=0:1000
  x(i+1)=feval(fun,(i/1000.0),mu);
end
plot(0:1/1000:1,x);
hold on;
plot([0 1],[0 1]);
xlabel('x(k)');
ylabel('x(k+1)');


y=zeros(1,nstep+2);
y(1)=xzero;

for i=1:nstep+1
  y(i+1)= feval( fun,y(i),mu);
end;

pmap = [];
y=[y(1:length(y)-1);y(2:length(y))];
for step=1:length(y)
  pmap=[pmap y(:,step) [y(2,step); y(2,step)]];
end;



pmap=pmap(:,1:length(pmap)-1);
plot(pmap(1,:),pmap(2,:));
plot(y(1,1),y(2,1),'*');


title('The return map:  * start');

hold off;