function out=simpson(f,a,b,N,varargin)
% SIMPSON('fun',A,B,N) numerically integrates the function
% 'fun' over the interval [A,B] using a fixed number of
% (2*N+1) function evaluations.
%
% For better, adaptive implemetations use QUAD or QUAD8.
% The advantage of SIMPSON is that is does NOT get HUNG UP
% at discontinuities etc ... it is plain dumb.
%
% All rights reserved: Matthias Kawski, September 30, 2000. 
% Contact: kawski@asu.edu, http://math.la.asu.edu/~kawski
%
dx=(b-a)/N;
xx=a+dx*[1:N-1];
out=sum([feval(f,[a,b],varargin{:}),...
      2*feval(f,xx,varargin{:}),...
      4*feval(f,[a,xx]+dx/2,varargin{:})])*dx/6;