% integrate.m
% 
% AREA = INTEGRATE(F,X0,X,H)
% 
%     INPUTS: F    --> Function to Integrate
%             X0   --> Lower Limit
%             X    --> Upper Limit
%             H    --> Accuracy of Integration (smaller is more accurate)
%
%    OUTPUTS: AREA --> Area under the curve bounded by the x-axis, the
%                      function, and the limits of integration.
% 
% This MATLAB function integrates any given algebraic function over a
% continuous and differentiable specified range to an accuracy inversely
% proportional to H.
% 
% While it is possible to have this function integrate over a discontinuous
% range, one must ensure that none of the numerical evaluation points
% coincide with a vertical asymptote.  To avoid this and get an
% approximation of an integral over an asymptote, add any random string of
% digits to the lower limit whose real-world value is less than H.
% 
%  BEHAVIOUR:      --> If H is less than 1e-6, it will default to 1e-6.
%                      This is to prevent ridiculously time-consuming
%                      numerical integrations.

% Copyright 06-17-2006 Brendan C. Wood -- Synchroverge

function area = integrate(f,x0,x,h)

def_h = 1e-6;

if x0 > x
    temp = x;
    x = x0;
    x0 = temp;
end

if h < 0
    fprintf('\nERROR: h cannot be less than zero!\n');
    return
else
    if h < def_h
        h = def_h;
    end
end

num_space = fix( ( x - x0 ) / h );

X = linspace(x0,(x-h),num_space);

Y = f(X);

A = Y * h;

area = sum(A);