Time for action - vectorizing the Monte Carlo integrator

Using the method above, the vectorized version of mcintgr is:

Code Listing 5.8
function I = mcintgrv(fun, a, b, mcloops) #1
#2
# Find maximum value of f #3
x = linspace(a,b); #4
f = feval(fun,x); #5
maxf = max(f); #6
#7
# Generating random arrays #8
r1 = rand(mcloops,1); #9
r2 = rand(mcloops,1); #10
#11
# Get random x and y values #12
l = b - a; #13
x = a + l.*r1; #14
y = maxf.*r2; #15
fx = feval(fun, x); #16
#17
# Counts the number of points that lie under the graph #18
counter = length(find(y<fx)); #19
#20
# The integral #21
I = counter/mcloops*maxf*l; #22
#23
endfunction #24

What just happened?

The vectorized function is called with the exact same inputs and outputs as mcintgr. The code should be straightforward to understand, perhaps except line 19. Here we first perform a comparison between the value of the function and the random points. This comparison operation results in a Boolean array which is passed to find, that returns...