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

x1=(a+h):h:b;
x2=(c+h):h:d;
[X1,X2]=meshgrid(x1,x2);
fx=feval(f, X1, X2);
q = sum(sum(fx*h*h));

