integrate(func, lb, ub, intfunction)

function result = integrate(func, lb, ub, intfunction)
    %%%%  Ariel Shwayder's integrate function %%%%%
    %%%% (c) 2008
    %%%%
    %%%% This is basically just a wrapper to make life easier in dealing
    %%%% with infinite bounds.
    %%%%
    %%%% In order to deal with infinite bounds it assumes that
    %%%%    (1/x^2)*f(1/x) is integratable, at least in a specific region.
    %%%%
    %%%% USAGE:
    %%%% result = integrate(func, lb, ub, [intfunction])
    %%%%    intfunction is optional, if none is given 'quad' will be used.
    
    if (nargin ~= 3 && nargin ~=4); error('Wrong # of arguments to integrate'); end;

    if nargin==3
        intfunction = @quad;
    end
    
    txffunc = @(x)(1./(x.^2)).*func(1./x);
    
    
    if (lb>ub)
        % Bounds are flipped
        result = -integrate(func,ub,lb,intfunction);
    elseif (lb==ub)
        result = 0;
    elseif (isfinite(lb) && isfinite(ub))
        % if both bounds not infinite, then can just do normal integral
        result = intfunction(func,lb,ub);
    elseif (lb==-inf && ub==inf)
        % We're trying to do -inf to inf
        result = intfunction(txffunc,-1,0)+intfunction(func,-1,1) ...
                    + intfunction(txffunc,0,1);
    elseif (isfinite(lb) && ub==inf)
        % We're trying to do a to inf
        if (lb>0)
            result = intfunction(txffunc,0,1/lb);
        else
            result = intfunction(func,lb,1)+intfunction(txffunc,0,1);
        end
    elseif (lb==-inf && isfinite(ub))
        % We're trying to do -inf to a
        if (ub<0)
            result = intfunction(txffunc,1/ub,0);
        else
            result = intfunction(txffunc,-1,0)+intfunction(func,-1,ub);
        end
    end