% bifurcation diagram for dy/dt=y(1-y)-h
clf
h=linspace(0,0.25);
y1=(1+sqrt(1-4*h))/2;
y2=(1-sqrt(1-4*h))/2;
%plot(h,y1,h,y2)
axis([0,.5,-0.2,1.2])
xlabel('h','FontSize',16)
ylabel('y','FontSize',16)
hold on
for hh=0:.05:.5
    plot([hh,hh],[-0.2,1.2],'k')
    if hh<=0.25
        yy1=(1+sqrt(1-4*hh))/2;
        yy2=(1-sqrt(1-4*hh))/2;
        plot([hh,hh],[yy1,yy2],'ok','MarkerSize',16,'MarkerFaceColor','k')
        plot([hh,hh],[-0.1,1.1],'vg','MarkerSize',16,'MarkerFaceColor','g')
        if hh<.25
            plot([hh],[0.5],'^g','MarkerSize',16,'MarkerFaceColor','g')
        end
    else
    plot([hh,hh,hh],[-0.1,0.5,1.1],'vg','MarkerSize',16,'MarkerFaceColor','g')
    end
    pause 
end
pause
plot(h,y1,h,y2)