parpool('local');

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

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

disp(res(end))
toc