N = 100000;


X = rand(N,1);
Y = rand(N,1);

R = sqrt(X.^2 + Y.^2);

ind = find(R < 1);

M = length(ind);

abs(M/N-pi/4)
figure(1)
plot(X,Y,'.k')
hold on
plot(X(ind),Y(ind),'.r')
hold off
axis equal
axis([0 1 0 1])