function bifurf(ra,rb,xa,xb,maxiter)
% Function evaluation is given by
%    bifurf(1,4,0,1,200)
%    bifurf(3.57,3.58,0.33,0.4,500000)
%
% modification of code in Matlab Guide, Higham and Higham, SIAM2005 
r=linspace(ra,rb,1001)';
ic=0.2:0.3:1;
pts=50;
rvals=repmat(r,length(ic),1);
xresults=zeros(length(rvals),pts);
x=kron(ic',ones(size(r)));

for k=2:maxiter
    x=rvals.*x.*(1-x);
end
for k=1:pts
   x=rvals.*x.*(1-x);
   xresults(:,k)=x;
end 

plot(rvals,xresults,'b.','MarkerSize',1)
axis([ra,rb,xa,xb]);
xlabel('r');
ylabel('x');
title('Bifurcation for Logistic Map')