function y1 = fpu( t, y )
%
%   Usage:
%   y' = fpu( t, y )  or  H = fpu( y )  or  y0 = fpu()
%
%   Fermi-Pasta-Ulam problem, y = [q,p].
%   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.
%
m   = 7;
m2  = 2*m;
w   = [10 10 10 10000 10 10 10]';
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 % f(y)
   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 = [p -fq];
elseif nargin==0  % y0
    q  = [0:m2-1]/(2*(m2-1)); 
    p  = zeros(1,m2);  
    y1 = [q p];    
else
    error
end
return