area = gquad (fun,xlow,xhigh,mparts,bp,wf)
or
area = gquad (fun,xlow,xhigh,mparts,nquad)
or
area = gquad (fun,xlow,xhigh,mparts,bp,wf,y)
This function evaluates the integral of an externally
defined function fun(x) between limits xlow and xhigh. The
numerical integration is performed using a composite Gauss
integration rule. The whole interval is divided into mparts
subintervals and the integration over each subinterval
is done with an nquad point Gauss formula which involves base
points bp and weight factors wf. The normalized interval
of integration for the bp and wf constants is -1 to +1. The
algorithm is described by the summation relation
x=b j=n k=m
integral( f(x)*dx ) = d1*sum sum( wf(j)*fun(a1+d*k+d1*bp(j)) )
x=a j=1 k=1
where bp are base points, wf are weight factors
m = mparts, and n = length(bp) and
d = (b-a)/m, d1 = d/2, a1 = a-d1
The base points and weight factors must first be generated
by a call to grule of the form [bp,wf] = grule(nquad)
Optional argument, nquad, is used if the Gauss points and weights
have not been previously calculated.
Optional argument, y, is used if the function, fun is a function
of x and y. fun(x,y) will be integrated over the range in x for
the constant, y.
