Matlab Calculo de una integral por números aleatorios

[isacampo]

AL.6. Los números aleatorios pueden emplearse para calcular integrales definidas y en particular para calcular áreas y volúmenes. El procedimiento es parecido al que hemos visto antes para obtener el número π.
Para el caso del área entre una curva y el eje de coordenadas,
Generamos un punto x aleatorio en el intervalo de integración [a,b].
Generamos un punto y aleatorio en el intervalo [0,ymax]. Donde ymax, representa el valor máximo que toma la función que queremos integrar en el intervalo [a,b].
Comprobamos el valor que toma la función en el punto x generado. Si el valor resultante es mayor el valor de y obtenido aleatoriamente, incrementamos un contador en una unidad. Se trata de puntos que están por debajo de la curva y, por tanto caen dentro del área comprendida entre la curva, los límites de integración y el eje x.
Repetimos este proceso un número grande de veces.
Al final, el número de puntos obtenidos por debajo de la curva frente al número total de puntos generados está en la misma proporción que el área encerrada por la curva frente al área del rectángulo formado por los límites de integración e ymax,


Escribe un programa que calcule empleando números aleatorios, la integral

Comprueba que la precisión obtenida mejorar con el número de puntos aleatorios generados para obtener la integral.

Creé este programa, pero no me da el resultado esperado. ¿Alguien puede ayudarme?

function[]=calcula_integral(n)
x=1+(2-1).*rand(1,n);
y=0+(1-0).*rand(1,n); %la y max para la funcion dada es 1.
l=length(x);
f=[];

for i=1:l
fu(i)=1/i;
f(i)=1./y(i);
end

for j=1:l
if fu(j)>f(j)
nfuera(j)=f(j);
elseif fu(j)<f(y)
ndentro(j)=f(j);
end
end

nf=0;
for k=1:l
if nfuera(k)~=0;
nf=nf+1;
end
end
nd=0;
for k=1:l
if ndentro(k)~=0;
nd=nd+1;
end
end

integral=1*(2-1)*(nd/n)
end

El ejercicio de calcular el número pi del que habla el ejercicio es el siguiente:

function[]=genera_pi(n)
x=rand(1,n);
y=rand(1,n);
l=length(x);
for i=1:l
d(i)=sqrt(x(i).^2+y(i).^2); %modulo de las coordenadas del punto
end

for j=1:l
if d(j)>1
nfuera(j)=d(j);
else
nfuera(j)=0;
end
end
for j=1:l
if d(j)<=1
ndentro(j)=d(j);
else
ndentro(j)=0;
end
end
%quita los ceros que hay en los vectores creados
nf=0;
for k=1:l
if nfuera(k)~=0;
nf=nf+1;
end
end

nd=0;
for k=1:l
if ndentro(k)~=0;
nd=nd+1;
end
end

numero_pi=4*(nd/n);
fprintf(1,'El valor de pi es %f \n',numero_pi)
end