function [Is,It]=integrali2d(f,a,b,c,d,h,k)
%funkcija nalazi vrednost integrala funki=cije f(x,y)
%na pravougaoniku [a,b] times [c,d] koristeci korak h na x osi,
% tj korak k na y osi
% test:
%[Is,It]=integrali2d(@(x,y) sin(x.^2+y.^2),0,0.6,0,0.4,0.3,0.2)
%Ugradjena funkcija: integral2(f,xmin,xmax,ymin,ymax), gde je f anonimna
%funkcija dve promenljive
%Za nas zadatak: integral2(@(x,y) sin(x.^2)+sin(y.^2),0,0.6,0,0.4)

xx=a:h:b;
n=length(xx);
yy=c:k:d;
m=length(yy);

F=ones(m,n);
for i=1:m
    for j=1:n
        F(i,j)=f(xx(j),yy(i));
    end
end


%%%trapezna formula
At=ones(m,n);
At(1,2:n-1)=2;
At(m,2:n-1)=2;
At(2:m-1,1)=2;
At(2:m-1,n)=2;
At(2:m-1,2:n-1)=4;

At
It=h*k*sum(sum((At.*F)))/4;

%%Simpsonova formula
if mod(n,2)~=0 && mod(m,2)~=0
    As=ones(m,n);
    
    pom=4*ones(1,m-2); 
    pom(2:2:end)=2;%pomocni vektor za prvu i poslednju kolonu
    As(2:m-1,1)=pom;
    As(2:m-1,n)=pom;

    pom2=4*ones(1,n-2); 
    pom2(2:2:end)=2;%pomocni vektor za prvu i poslednju vrstu
    As(1,2:n-1)=pom2;
    As(m,2:n-1)=pom2;
  
    As(2:2:m-1,2:2:end-1)=16;
    As(3:2:m-1,2:2:end-1)=8;
    As(2:2:m-1,3:2:end-1)=8;
    As(3:2:m-1,3:2:end-1)=4;
    
    As

    Is=h*k*sum(sum(As.*F))/9;
else
    disp('Broj cvorova na obe ose mora biti neparan!')
end