function y1 = argon2( t, y )
%
%   Usage:
%   q'' = F*argon(t,q)  or  H = argon([q,p])  or  F = argon() == 1/M;
%
%   Frozen argon crystal problem in second order form.
%
%   See pages 144-145 of
%   L.Brugnano, F.Iavernaro.
%   Line Integral Methods for Conservative Problems.
%   Chapman and Hall/CRC - 2015
%   Series: Monographs and Research Notes in Mathematics.
%   
kB = 1.380658e-23;   % Boltzmann constant
ee = 119.8*kB;
N  = 7;              % number of atoms
M  = 66.34e-27;      % atom  mass
ss = 0.341;
s2 = ss^2;
if nargin==0         % F
    y1 = 1/M;        % it is a scalar matrix
    
elseif nargin==1  % Hamiltomian
   y  = t; 
   [m,n] = size(y);
   q  = y(:,1:n/2);
   p  = y(:,n/2+1:n);
   q1 = q(:,1:2:end);
   q2 = q(:,2:2:end);
   p1 = p(:,1:2:end);
   p2 = p(:,2:2:end);
   Hp = (1/(2*M))*sum( ( p1.^2+p2.^2 ).' ).';
   y1 = 0;
   for i = 2:n/4
       for j=1:i-1
           rij = s2./( ( q1(:,i)-q1(:,j) ).^2 + ( q2(:,i)-q2(:,j) ).^2 );
           y1 = y1 + ( rij.^6 - rij.^3 ); 
       end
   end    
   y1 = Hp + (4*ee)*y1;
   
else             % function
   [m,n] = size(y); n = 2*n;
   q1    = y(:,1:2:end);
   q2    = y(:,2:2:end);
   Hq    = zeros(m,n/2);
   ind   = [1 2];
   for i = 1:n/4
       for j = [1:i-1 i+1:n/4]
           rij = s2./( ( q1(:,i)-q1(:,j) ).^2 + ( q2(:,i)-q2(:,j) ).^2 );
           fij = ( (6/s2)*rij.^7 - (3/s2)*rij.^4 );
           Hq(:,ind) = Hq(:,ind) -[fij.*( q1(:,i)-q1(:,j) ) fij.*( q2(:,i)-q2(:,j) )];
       end
       ind = ind +2;
   end    
   y1 = (-2*4*ee)*Hq;    
end    
