function q = my_quad2D(f, a, b, c, d, h)

Ni = ceil((b-a)/h);
Nj = ceil((d-c)/h);
q = 0;

for i = 1:Ni 				% loop över samtliga delintervall - x-led
  for j = 1:Nj 				% loop över samtliga delintervall - y-led
    x1i = a + i*h; 			% delkvadratens högra ändpunkt
    x2j = c + j*h; 			% delkvadratens övre ändpunkt
    q = q + feval(f, x1i, x2j)*h*h; 	% ``right-hand rectangle rule''
  end
end