tic
broj_tacaka=1000;
delta=1; % parametar za sirinu uniformne raspodele
tacke=zeros(2,broj_tacaka);
brojac=0;
tacke(:,1)=[0;0];
stara_tacka=tacke(:,1);
for i=2:broj_tacaka
 nova_tacka=stara_tacka+delta*(rand(2,1)-0.5);
stara_p=exp(-(stara_tacka(1)^2)*(stara_tacka(2)^2)-stara_tacka(1)^2-stara_tacka(2)^2);
nova_p=exp(-(nova_tacka(1)^2)*(nova_tacka(2)^2)-nova_tacka(1)^2-nova_tacka(2)^2);
 if rand<(nova_p/stara_p)
  stara_tacka=nova_tacka;
  brojac=brojac+1;
 end
 tacke(:,i)=stara_tacka;   
end
uzorak=tacke(1,:).^2;
integral=sum(uzorak)/broj_tacaka;
figure
plot(tacke(1,:),tacke(2,:),'-o')
xlabel('x');
ylabel('y');
a=integral2(@(x, y) x.*x.*exp(-x.*x-y.*y-x.*x.*y.*y), -20, 20, -20, 20);
b=integral2(@(x, y) exp(-x.*x-y.*y-x.*x.*y.*y), -20, 20, -20, 20);
fprintf ( 1, 'Vrednost integrala dobijena numerickom integracijom na mrezi je: %f \n', a/b);
fprintf ( 1, 'Vrednost integrala dobijena Metropolisovom MC metodom je: %f \n', integral);
fprintf ( 1, 'Procenat prihvacenih tacaka: %f %% \n', brojac/(broj_tacaka-1)*100);
toc