function y1 = henon2(t,q)
%
%   Usage:
%   q'' = M*henon2( t, q )   or   M = henon2()   or  H = henon2([q,p]) ==1
%
%    Henon-Heiles problem in second-order form. 
%
%   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     % M
   y1 = 1;  % when it is a scalar matrix, the scalar suffices      
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   
   y1 = [-q(:,1).*(1+2*q(:,2)) -q(:,2)-q(:,1).^2+q(:,2).^2];
end
return