function bifurcation

% This program takes a one-dimensional system with a parameter r,
% plots the bifurcation diagram over a specified range,
% It also plots the phase portrait, and allows the user to
% interactively change the value of r to see how the bifurcation
%takes place

%NOTES ON USAGE:
%(1)  the bounds for the phase potrait axis and parameter must be
%     changed for each problem
%(2)  one MUST define f(x,r) using ARRAY operations 
%     .*  ./  .^ instead of *, /, ^.  For example,
%  f(x,r) = rx - x^2/(1+x)  would be entered as
%
%  f = @(x,r)(r*x - (x.^2)./(1+x)); 
%
%  If you want to understand more about this, use the help feature.
% (3)  To exit the interactive program, hit Control+C
%    (Sorry, some day I'll get around to putting in a quit button)


% The following are problem-specific; CHANGE THESE

axisbounds = [0 1  -.2  .5];   %bounds for phase portrait 

rbounds = [ 0 1];   %minimum and maximum for parameter r

%%%%%  This is the fishery model problem (Strogatz 3.7.4) with a = .5

a = .5;  % this parameter is fixed for purposes of investigating the bifurcation in h

f = @(x,h)(x.*(1- x) - h*x./(a + x)); %DONT FORGET to use .*  ./  .^ instead!


%%%%%%% THE rest does not need to be edited
%%%%%%%  plot bifurcation diagram
figure(1); clf;
r = (rbounds(1)+rbounds(2))/2;  % initial value for r
%plot bifurcation diagram
figure(1);  subplot(2,1,1);
for rr = [rbounds(1):0.01:rbounds(2)];
 x = [axisbounds(1):0.01:axisbounds(2)];
 fx = f(x,rr);
 [fz,sta,nz]  = findzeros(x,fx);  
 figure(1);  
 title('bifurcation diagram','fontsize',16); 
 xlabel('parameter', 'fontsize',14); ylabel('fixed points', 'fontsize',12);
 hold on;
 for j=1:nz,
    if (sta(j)>0) plot(rr,fz(j),'.'); else plot(rr,fz(j),'o');   end;      
 end;
 axis([rbounds axisbounds(1:2)])
end;

% now loop for interactive phase portraits
h = line([r r],axisbounds(1:2),'linestyle','--');
while(1)
 x = [axisbounds(1):0.01:axisbounds(2)];
 fx = f(x,r);
 figure(1);  subplot(2,1,2);
 plot(x,fx);  xlabel('x'); ylabel('xdot'); hold on;
 [fz,sta,nz]  = findzeros(x,fx);  %find zeros and stability of f;
 for j=1:nz,
     if (sta(j)>0) plot(fz(j),0,'.'); else plot(fz(j),0,'o');  end;  
 end;
 axis(axisbounds);
 title( ['Interactive phase portrait:  r = ' num2str(r)],'fontsize',20);
 line( [axisbounds(1) axisbounds(2)],[0 0],'color','black');
 hold off;
 text(axisbounds(1)+.05,.05,'click to decrease r');
 text(axisbounds(2),.05,'click to increase r','HorizontalAlignment','right');
 pause(.1);
 [X Y] = ginput(1);  %get mouse position
 if ( X < (axisbounds(1)+axisbounds(2))/2)  r = r - .05; else r =r +.05; end;
 delete(h);  subplot(2,1,1);
 h = line([r r],axisbounds(1:2),'linestyle','--');
end;

%-----------------------------------------------------------
%function to find zeros, stability of f(x,r) 
function [fzeros,sta,nz] = findzeros(x,fx);
   df = -.5*(sign(fx(2:end)).*sign(fx(1:end-1))-1).*sign(fx(2:end)-fx(1:end-1));
   nz = nnz(df);
   i = find(df);
   sta = df(i);
   fzeros = x(i);
%-------------------------------------------------------------