clear; format long g m = 1e7; A = rand(m,2); f = ( A(:,2) <= sin(A(:,1)./A(:,1)) ); n = sum(f); p = n / m % p = 0.8413388 quadl('sin(x)./x',0,1) % ans = 0.946083070367176