function I = e4_Simpson(f,a,b,n,varargin)
% Simpson: composite 1/3 Simpson rule 
% I = Simpson(func,a,b,n,pl,p2,...):
% 
% input:
% f = name of function to be integrated
% a, b = integration limits
% n = number of segments (default = 100)
% pl,p2,... = additional parameters used by func
% output:
% I = integral estimate

if floor(n/2)*2 ~= n, error('n should be even'), end

if nargin<3,error('at least 3 input arguments required'),end
if ~(b>a),error('upper bound must be greater than lower'),end
if nargin<4|isempty(n),n=100;end
h = (b - a)/n;
s1 = 0;
for i = 1:2:n-1
x = a +i*h;
s1 = s1 + 4*f(x,varargin{:});
end
s2 = 0;
for i = 2:2:n-2
x = a +i*h;
s2 = s2 + 2*f(x,varargin{:});
end
s =f(a,varargin{:}) + s1 + s2 + f(b,varargin{:});
I = h/3 * s;