function [I] = Int2D(fxy,a,b,c,d,h)
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
f=inline(fxy);
isx=(b-a)/(2*h);
isy=(d-c)/(2*h);
A=zeros(isx,isy);
for m=1:isy
for n=1:isx
A(n,m)=h^2/9*(f(a,c)+f(a+2*h,c)+f(a,c+2*h)+f(a+2*h,c+2*h)+4*(f(a+h,c)+f(a,c+h)+f(a+2*h,c+h)+f(a+h,c+2*h))+16*f(a+h,c+h));
a=a+2*h;
end
a=a-2*h*isx;
c=c+2*h;
end
l=0;
for m=1:isy
for n=1:isx
l=l+A(n,m)
end
end
l
end
