rombergQuadrature

% [x, err] = rombergQuadrature(fun,tSpan,tol)
%
% Compute the integral(fun), over the domain tSpan, to an accuracy of tol,
% using Romberg quadrature. Fully vectorized.
%
% Good for high-accuracy quadrature over smooth vector functions.
%
% If necessary, this function will automatically sub-divide the interval to
% achieve the desired accuracy. This should only occur when fun is stiff or
% non-smooth.
%
% INPUTS:
% fun = vector function to be integrated
% dx = fun(t)
% t = [1, nt] = time vector
% dx = [nx, nt] = function value at each point in t
% tSpan = [tLow, tUpp] = time span (domain) for integration
% tol = [nx,1] = desired error tolerance along each dimension
%
% OUTPUT:
% x = [nx,1] = integral along each dimension
% err = [nx, 1] = error estimate along each dimension
%
% NOTES:
% algorithm from:
% http://www.math.usm.edu/lambers/mat460/fall09/lecture29.pdf
%