function I=trapez2d(fun,a,b,c,d,n)
%I=trapez2d(fun,a,b,c,d,n)
%izracuna dvojni integral
%int_a^b dx int_c^d dy f(x,y)

%computes the double integral
%int_a^b dx int_c^d dy f(x,y)
%using the trapese rule

x=linspace(a,b,n);%vozli v x smeri // knots in the x direction
y=linspace(c,d,n);%vozli v y smeri // knots in the y direction
hx=x(2)-x(1);%dolzina podintervala v x smeri // length of subinterval in x direction
hy=y(2)-y(1);%dolzina podintervala v y smeri // length of subinterval in y direction
[X,Y]=meshgrid(x,y);%matriki vrednosti za x in y na mrezi tock // meshgrid of points
fval=feval(fun,X,Y);%vrednosti funkcije na mrezi tock // compute function values on the meshgrid
w=[1 2*ones(1,n-2) 1];%utezi za 1-dimenzionalen primer// weights for 1 dimensional integration method
w=w'*w;%matrika utezi za 2d je tenzorski produkt utezi za 1d // weight fro 2d is the tensor produkt of weights in 1d
I=sum(sum(w.*fval))*hx*hy/4;%priblizek za integral // approximation of the integral