Ok, now that we understand the question, it is easy. Just use trapz, TWICE.
To apply the trapezoidal rule, you will create a matrix of 4x4 elements. Then call trapz twice, once on each dimension of the array. I won't do your problem here, since this seems too much likely to be homework.
But for example, suppose I wanted to integrate the function
z = x.^3 + x.*y^2
over the domain (x,y) = [20,85] X [10,98], where we will sample the function at that set of nodes? This is VERY different from what you wanted to call a discontinuous point set of some sort.
Zfun = @(X,Y) X.^3 + X.*Y.^2;
[x,y] = meshgrid(Xnodes,Ynodes)
intest = trapz(Xnodes,trapz(Ynodes,z,1),2)
How accurate is that estimate? Remember, this is just trapezoidal integration, over a VERY coarse set of nodes.
int(int(X.^3 + X.*Y.^2,X,20,85),Y,10,98)
So the trapezoidal estimate was within 10% of the exact value. Probably as good as we could do on that set of nodes.