f = @(x,y) 2*x.^2 + y.^2;
integrate_in_2d( f, [0,2], [0,1] )
