function compute_bifn2
% Collection of commands to compute bifurcation curves
% This script computes the subcritical pitchfork bifurcation of func2

% First plot the vector field, to estimate fixed point locations
xx = (-1:0.01:1);
figure(1)
clf    % clear figure 1
% Try r = -0.2 (below bifurcation point)
plot(xx,func2(xx,-0.2),'b',xx,zeros(size(xx)),'--')
pause  % hit any key to continue

% There seem to be zeros near x = -0.5, x = 0 and x = 0.5
% Use fzero to find them accurately, bracketing the intervals
% There are also fixed points near x = -0.9 and x = 0.9, but we will find
% these by continuing the bifurcation curves from the other nonzero 
% equilibria.
% Note that func2(x,-0.2) explicitly specifies r = -0.2
x1 = fzero('func2(x,-0.2)',[-0.7 -0.2]);
x2 = fzero('func2(x,-0.2)',[-0.2 0.2]);
x3 = fzero('func2(x,-0.2)',[0.2 0.7]);
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
% ti and tf to choose.

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

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