%
% Simpson composite rule
%
% fun = funtion to integrate
% a   = left boundary
% b   = right boundary
% n   = number of interval (MUST be even!)
%
% end = is the last array index
%
% (h/3)*(f(a)+f(b))+ (h/3)*(4*f(x(2))+2*f(x(3))+...+4*f(x(N-1))
%
%
function RES = simpson( fun, a, b, n )
  h     = (b-a)/n ;
  nodes = [a:h:b] ;
  yvals = feval( fun, nodes ) ;
  RES   = yvals(1)+yvals(end) ; % f(a) + f(b)
  RES   = RES + 4*sum(yvals(2:2:end-1)) ;
  RES   = RES + 2*sum(yvals(3:2:end-2)) ;
  RES   = (h/3)*RES ;
end
