function compute_bifn
% Collection of commands to compute bifurcation curves
% This file computes the supercritical pitchfork bifurcation of func1

% First plot the vector field, to estimate fixed point locations
xx = (-5:0.1:5);
figure(1)
clf    % clear figure
% Try r = 2
plot(xx,func1(xx,2),'b',xx,zeros(size(xx)),'--')
pause  % hit any key to continue

% There seem to be zeros near x = -2, x = 0 and x = 2
% Use fzero to find them accurately, bracketing the intervals
% Note that func1(x,2) explicitly specifies r = 2
x1 = fzero('func1(x,2)',[-5 -0.5]);
x2 = fzero('func1(x,2)',[-0.5 0.5]);
x3 = fzero('func1(x,2)',[0.5 5]);
fprintf('x1 = %6.3f\n', x1)
fprintf('x2 = %6.3f\n', x2)
fprintf('x3 = %6.3f\n', x3)
pause

% Now use these values as initial values to trace out the bifurcation 
% curves depending on r.  We will integrate forwards and backwards, to 
% tf > 0 and ti < 0; it might take some experimentation to see what times
% to integrate to.

figure(2)
clf
hold on
bifurcationplot('vfield1',x1,2,-1,10)
bifurcationplot('vfield1',x2,2,-10,0.8)
bifurcationplot('vfield1',x3,2,-1,10)

% The curves beginning at r = 2 cannot get past the bifurcation point
% at r = 1; so explicitly compute the curve below the pitchfork bifurcation
% (which we know has a unique branch)
x4 = fzero('func1(x,0)',[-1 1]);
bifurcationplot('vfield1',x4,0,-10,1)