% [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
% 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
% x = [nx,1] = integral along each dimension
% err = [nx, 1] = error estimate along each dimension
% algorithm from:
Now the rombergQuadrature automatically detects a non-smooth integrand and sub-divides the interval to achieve the desired accuracy.