clear
N = 100;
Mrep = 1000;
pi_exp = zeros( Mrep, 1 );
for j=1:Mrep
    pi_exp(j) = 4*mean( sum( rand(N,2).^2,2 )<1 );
end
figure, plot( pi_exp, '*' )
std( pi_exp )
