function y1 = fpu2( t, y )
%
%   Usage:
%   q'' = M*fpu2( t, q )  or  H = fpu2( [q p] )  or  M = fpu1() == 1
% 
%   Fermi-Pasta-Ulam problem  in second order form. 
%   7 hard springs with constants w = [10 10 10 1e4 10 10 10]. 
%
%   See pages 136-137 of
%   L.Brugnano, F.Iavernaro.
%   Line Integral Methods for Conservative Problems.
%   Chapman and Hall/CRC - 2015
%   Series: Monographs and Research Notes in Mathematics.
%
w   = [10 10 10 10000 10 10 10]';
m   = length(w);
m2  = 2*m;
w2  = w.^2/2;
uno = ones(m+1,1);
un2 = ones(m2,1)/2;
if nargin==1 % Hamiltonian
   [n,m4] = size(t);
   qo     = t(1:n,1:2:m2);
   qe     = t(1:n,2:2:m2);
   p      = t(1:n,m2+1:end);
   O      = zeros(n,1);
   y1     = p.^2*un2 +(qe-qo).^2*(w2/2) + ([qo O]-[O qe]).^4*uno;
elseif nargin==2 % M * nabla U(q)
   n  = length(t);
   qo = y(1:n,1:2:m2);
   qe = y(1:n,2:2:m2);
   O  = zeros(n,1);
   un = ones(n,1);
   dq = 4*([qo O]-[O qe]).^3;
   %p  = y(1:n,m2+1:end);
   fq = zeros(n,m2);
   fq(:,1:2:m2) = (qo-qe).*(un*w2') +4*( qo -[O qe(:,1:end-1)] ).^3 ;
   fq(:,2:2:m2) = (qe-qo).*(un*w2') -4*( [qo(:,2:end) O] -qe ).^3 ;
   y1 = -fq; 
elseif nargin==0  % M 
    y1 = 1;  % scalar is also allowed when the matrix is scalar
else
    error
end
return