Numerical integration using adaptive Simpson's rule.

## Inputs

func
The function to integrate.
The input is a scalar and the output is a scalar.
A function handle or anonymous function is also supported.
a
Lower integration limit.
Type: double
Dimension: scalar | vector | matrix
b
Upper integration limit.
Type: double
Dimension: scalar | vector | matrix
abstol
Absolute tolerance (default: sqrt(eps) or about 1.0e-8).
Type: double
Dimension: scalar

## Outputs

area
The estimated area.
Dimension: scalar | matrix
count
Number of function evaluations.
Dimension: scalar | matrix

## Examples

Single interval:
``````function y = Integrand(x)
y = 13 * (x - x^2) * exp(-3*x/2);
end

[area,count] = quadv(@Integrand, 0, 4, 1.0e-5)``````
``````area = -1.54878858
count = 69``````

Multiple intervals with an anonymous function:

``[area,count] = quadv(@(x) sqrt(x), [0, 1], [1, 2])``
``````area = [Matrix] 1 x 2
0.66667  1.21895
count = [Matrix] 1 x 2
309  29``````

quadv recursively bisects each interval until the improvement from the bisections falls below the absolute tolerance. Intervals are assumed to have finite bounds.

The maximum number of function evaluations is 10,000, and the minimum interval is 1.0e-12.

To pass additional parameters to a function argument, use an anonymous function.

quadv does not yet support func returning a matrix to support multiple simultaneous integrands.