error: vol = gquad2d(fun,xlow,xhigh,ylow,yhigh,bpx,bpy,wfxy) or vol = gquad2d(fun,xlow,xhigh,ylow,yhigh,nquadx,nquady) This function evaluates the integral of an externally defined function fun(x,y) between limits xlow and xhigh and ylow and yhigh. The numerical integration is performed using a Gauss integration rule. The integration is done with an nquadx by nquady Gauss formula which involves base point matrices bpx and bpy and weight factor matrix wfxy. The normalized interval of integration for the bpx, bpy and wfxy constants is -1 to +1 (in x) and -1 to +1 (in y). The algorithm is described by the summation relation x=b j=nx k=ny integral( f(x)*dx ) = J*sum sum( wfxy(j,k)*fun( x(j), y(k) ) ) x=a j=1 k=1 where wfxy are weight factors, nx = nquadx = number of Gauss points in the x-direction, ny = nquady = number of Gauss points in the y-direction, x = (xhigh-xlow)/2 * bpx + (xhigh+xlow)/2 = mapping function in x, y = (yhigh-ylow)/2 * bpy + (yhigh+ylow)/2 = mapping function in y, and J = (xhigh-xlow)*(yhigh-ylow)/4 = Jacobian of the mapping. The base points and weight factors must first be generated by a call to grule of the form [bpx,bpy,wfxy] = grule2d(nquadx,nquady) The first form of gquad2d is faster when used several times, because the points and weights are only calculated once. The second form of gquad2d is usefull if it is only called once (or a few times).