function T = trapez(f,a,b,N)
% Composite trapezoidal 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,N+1); % Initialize gridpoints
fv = f(x); % row vector of function values
fv(1) = fv(1)/2; fv(end) = fv(end)/2; % special treatment of endpoints
T = (b-a)*sum(fv)/N; % summation