function z=integrateFunc(f,x1,x2,dx)

% Given a function f(x), returns its definite integral between x1 and x2. 
% Uses a linear interpolation approximation for y(x) between points.
% Input:
%   f(x)   : function to be integrated (not a vector of values)
%   x1, x2 : lower and upper bounds of integration
%   dx     : integration step
% Output:
%   z      : integral of f(x)
% Usage:
%   F = integrateFunc( @myFunc, 0, 2, 0.01 )
% J.M.Soler. Nov.2010

nx = ceil((x2-x1)/dx);      % number of integration intervals
dx = (x2-x1)/nx;            % so that last point is exactly x2
x = (x1 : dx : x2);         % integration points
y = f(x);                   % function value at integration points
nx = numel(x);              % number of integration points
z = (0.5*y(1) + sum(y(2:nx-1)) + 0.5*y(nx))*dx;   % trapezoidal rule
end % function
