function bifurcationplot(fieldname,x0,r0,ti,tf)
% Computes bifurcation curves near r0, x0

options = odeset('AbsTol',1.e-12,'RelTol',1.e-8);
init = [r0 x0]';

% First integrate forwards from 0 to tf > 0
tspan = [0 tf];
[T,Y] = ode45(fieldname,tspan,init,options);
R = Y(:,1); X = Y(:,2);
plot(R, X, 'r-')
xlabel('r')
ylabel('x')
title(['Bifurcation curve : ' , fieldname])
hold on

% Now integrate backwards from 0 to ti < 0
tspan = [0 ti];
[T,Y] = ode45(fieldname,tspan,init,options);
R = Y(:,1); X = Y(:,2);
plot(R, X, 'b--')

    