gquad

  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.