clear
clf
hold on
abeg=input(' Enter the value of a to begin: ');
astep=input(' Enter the step in a: ');
afin=input(' Enter the final a: ');
ntot=input(' Enter total no. of iterations per a value: ');
a=abeg:astep:afin;    %Vector of all desired a values
na=length(a);         %Number of desired a values
x(1,1:na)=0.501;      %Choose generic initial condition, same for all a
for n=1:ntot
    x(n+1,:)=a.*x(n,:).*(1-x(n,:)); %Iterate all the a values at once
end
for n=ntot-80:ntot+1  %Plot only last 82 points (eliminates unwanted transients)
  plot(a,x(n,:),'r.','MarkerSize',2)
end
xbif=1+sqrt(6); % Value of a where two-cycle loses stability
yplot=0.3:.01:1;                %
xplot=ones(size(yplot))*xbif;   % Plot dashed vertical line at xbif
plot(xplot,yplot,'b--')         %