To compute U(i,j) and L(i,j) for each pair i and j, you will presumably have to use matlab's 'fmincon' function in the Optimization Toolbox which finds the minimum (infimum) of a function, f, subject to constraints x(i-1)<=x<=x(i) and y(j-1)<=y<=y(j), using -f in the case of U(i,j) to get the maximum (supremum) of f. You can read the details about 'fmincon' at:
http://www.mathworks.com/help/optim/ug/fmincon.html
Otherwise, I see no further problems in implementing the necessary code. As you have indicated, the summation can be accomplished using the 'sum' function with U and L being generated as N-by-M arrays.
Note that if f is a continuous function, then as the subdivisions in x and y become finer, the value of the double summation must approach zero, since each difference U(i,j)-L(i,j) will be approaching zero.
[The statement that N and M are "real" should actually state that they are integers. Otherwise the summation would make no sense.]
