TrapezoidalSum2D.m

function Int = TrapezoidalSum2D(f,x,y)

    xa = x(1:end-1);
    xb = x(2:end);
    ya = y(1:end-1);
    yb = y(2:end);
    dx = ones(length(xa),1)*(xb - xa);
    dy = (yb - ya)'*ones(1,length(ya));
    darea = dx.*dy;
    fa = f(1:end-1,1:end-1);
    fb = f(1:end-1,2:end);
    fc = f(2:end,1:end-1);
    fd = f(2:end,2:end);    
    Int = sum( sum( darea.*(fa + fb + fc + fd)./4 ) );
end