quad2dg

error:  int = quad2dg('Fun',xlow,xhigh,ylow,yhigh)
or
        int = quad2dg('Fun',xlow,xhigh,ylow,yhigh,tol)

This function is similar to QUAD or QUAD8 for 2-dimensional integration,
but it uses a Gaussian quadrature integration scheme.  
     int     -- value of the integral
       Fun     -- Fun(x,y) (function to be integrated)
       xlow    -- lower x limit of integration
       xhigh   -- upper x limit of integration
       ylow    -- lower y limit of integration
       yhigh   -- upper y limit of integration
       tol     -- tolerance parameter (optional)
Note that if there are discontinuities the region of integration 
should be broken up into separate pieces.  And if there are singularities,
a more appropriate integration quadrature should be used 
(such as the Gauss-Chebyshev for a specific type of singularity).