function [ xd ] = fsim( x )

%  evaluate x dot = f(x)
%  This one is for a one dimensional string of n springs
%  x = ( q_1,...,q_n,p_1,...,p_n), position and momentum variables
   
   n = length(x)/2;

   xd = zeros(1,2*n);
   
   for i = 1:n
       xd(i) = x(i+n);   %  q_i dot = p_i
   end
   
   for i = 2:(n-1)
       sl      = x( i ) - x(i-1);   %  string stretch to the left  of particle i
       sr      = x(i+1) - x( i );   %  string stretch to the right of particle i
       xd(i+n) = sr - sl;           %  spring constant = 1, force = stretch.
   end
   
   sl = x(1);          %  string stretch to the left  of particle i, assume x(0) = 0
   sr = x(2) - x(1);   %  string stretch to the right of particle i
   xd(1+n) = sr - sl;  %  spring constant = 1, force = stretch.
   
   sl      = x( n ) - x(n-1);   %  string stretch to the left  of particle n
   sr      =        - x( n );   %  string stretch to the right of particle n, x(n+1) = 0
   xd(n+n) = sr - sl;           %  spring constant = 1, force = stretch.
   
end
