function subpf(pffun,x0,r0,rf,nu)
 global epsilon ;   % a parameter for the slow evolution of r in pffun

 if (nargin < 4)
     epsilon = 1.0e-4;
 else
     epsilon = nu;
 end
  
 options = odeset('AbsTol',1.e-12,'RelTol',1.e-6,'OutputFcn','odephas2');
 % The option 'OutputFcn','odephas2' plots the first two components of the
 % solution after every timestep.
 
 init=[r0 x0]';
 t0 = 0;
 tf = 2*pi*(rf-r0)./epsilon;
 tspan=[t0 tf];
 
 figure(1);
 [T,Y]=ode15s(pffun,tspan,init,options);
 
 R = Y(:,1);  X = Y(:,2);
 
 figure(2)
 clf
 plot(R,X,'b')
 axis([min(R) max(R) -max(abs(X)) max(abs(X))])
 
 hold on
 xlabel('r')
 ylabel('x')
 title(['Subcritical pitchfork bifurcation : ' , pffun])
 
% Plot the bifurcation curves as well
  if strcmp(pffun,'mysubpf')
    plot([min(R) max(R)],[0 0],'r')
    plot(X.^4-X.^2,X,'r')
    plot(X.^4-X.^2,-X,'r')
  end

  return
% Comment out return to plot x and r vs t
    figure(3)
    clf
    hold on
    plot(T,Y(:,2),'r')
    plot(T,Y(:,1),'b')
    xlabel('Time')
    ylabel('r and x')

