tic
N_int=100;
N_int_0=1e2;
N_int_f=1e7;

passo=ceil((N_int_f-N_int_0)/N_int);
res=zeros(N_int,1);
i=1;
for N=N_int_0:passo:N_int_f
    Z=rand(2,N);
    res(i)=sum((Z(1,:)-0.5).^2+(Z(2,:)-0.5).^2<0.5^2)/(N*0.5^2);
    i=i+1;
end

disp(res(end))
toc