function T = simpson(f,a,b,N)
% Composite Simpson rule on an equidistant grid for the numerivcal
% integration of a continuous function.
% a,b are the interval bounds, f must be a handle of type @(x) to a
% function f(x), N is the number of grid intervals.
% Note: the function f must support vectorization
x = linspace(a,b,2*N+1); % Initialize gridpoints
fv = f(x); % row vector of function values
fv(2:2:end-1) = 4*fv(2:2:end-1); % midpoints of grid intervals
fv(3:2:end-2) = 2*fv(3:3:end-2); % grid points
T = (b-a)*sum(fv)/N/6; % summation