Integration in 2D by the midpoint rule

I'm trying to integrate a function ,log(1+x+y)*beta(x,2.5,3.0)*beta(y,2,4.5) ,by the midpoint rule, both x and y belongs to [0,1]. However, the approximation is far off from the actual value which is 0.0006469 . Could anyone help me to find where the errors might be? Thank you so much!

if true
  % code
end
clear
syms x;
syms y;
%creating the subintervals over the interval 0:1 by 10
a = 0 ;
b = 1 ;
n = 10;
h = (b-a)/n;
%calculating the midpoint for each subinterval
y = linspace(0,1,n+1);
for i = 1:n
  midpoint(i) = y(i)+ h/2;
end
%First, integrating w.r.t. y over the range from 0 to 1 by the midpint rule
g =@(x) 0;
for i = 1:n
  g = @(x) g(x) + h * log(1+midpoint(i)+x)* betapdf(midpoint(i),2,4.5);
end

if true
  % code
end
%Integration w.r.t. x over the range from 0 to 1 by the midpoint rule
Int = 0;
x=midpoint;
for i = 1:n
  Int = Int + g(x(i))*betapdf(x(i),2.5,3);
end;
Int

Answer 1

[Torsten]

fun=@(x,y)log(1+x+y)*beta(x,2.5,3.0)*beta(y,2,4.5);
n = 10;
h = 1/n;
x = linspace(0,1,n+1);
y = linspace(0,1,n+1);
int=0.0;
for ii=1:n
for jj=1:n
int=int+h^2*fun(x(ii)+h/2,y(jj)+h/2);
end
end

Best wishes

Torsten.