function [ t, y, te, ye, E ] = Henon( y0, tf )
%HENON  Integration of the Henon-Heiles model locating q[1] = 0
%   solve the model on 0 \le t \le tf, locating all places q[1]=0 

    function dy = Hf( t, y )
        dy = zeros(4,1);
        p = y(1:2);
        q = y(3:4);
        dy(1:2) = -[q(1) + 2*q(1)*q(2); q(2) + q(1)^2-q(2)^2 ];
        dy(3:4) = [p(1); p(2)];
    end

    function [z, isterm, dir] = Poincare(t,y)
        z = y(3);
        isterm = 0;
        dir = 1;
    end
    % Tell solver about events function, and make tolerances TIGHT.
    opts = odeset('Events', @Poincare, ...
                  'Reltol', 100*eps, 'Abstol', 1.0e-14 );

    [ t, y, te, ye ] = ode113( @Hf, [0,tf], y0, opts );
    E = (y(:,1).^2+y(:,2).^2+y(:,3).^2+y(:,4).^2)/2 ...
        + y(:,3).^2.*y(:,4)-y(:,4).^3/3;
    
end

