function y1 = henon(t,y)
%
%   Usage:
%   y0 = henon    or   y' = henon(t,y)    or   H = henon(y)
%
%   y = [q1 q2 p1 p2]. Henon-Heiles problem, first order formulation. 
%
%   See pages 137-139  of
%   L.Brugnano, F.Iavernaro.
%   Line Integral Methods for Conservative Problems.
%   Chapman and Hall/CRC - 2015
%   Series: Monographs and Research Notes in Mathematics.
%
if nargin==0     % initial condition
   y1 = [0 -1/8 99/1e4 -sqrt( (3/4)^4-1e-4 )];   
elseif nargin==1 % Hamiltonian
   q  = t(:,1:2);
   p  = t(:,3:4);
   y1 = ( p(:,1).^2+p(:,2).^2+q(:,1).^2+q(:,2).^2 )/2 + ...
        q(:,1).^2.*q(:,2)-q(:,2).^3/3;
else             % vector field 
   q = y(:,1:2);
   p = y(:,3:4);
   y1 = [p -q(:,1).*(1+2*q(:,2)) -q(:,2)-q(:,1).^2+q(:,2).^2];
end
return