As user1362215 points out, your function should be contained in the rectangle. You get closer to the solution if you increase n. Here is a vectorised solution. Results are in the range.
# Hit and miss
f <- function(x) x ^ (-0.5)
n <- 1000000
a <- 0.01
b <- 1
#ceiling(max(f((seq(0.01,1,by=0.001)))))
#[1] 10
set.seed(5)
x <- runif(n,a,b)
y <- 10*runif(n,0,1)
R <- sum(y < f(x))/n
(b-a)*10*R
#[1] 1.805701
# Repeat a few times to look at the distribution
set.seed(5)
n <- 100000
r <- replicate(1000,sum(10*runif(n,0,1) < f(runif(n,a,b)))/n *(b-a)*10)
hist(r)
summary(r)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1.755 1.792 1.800 1.800 1.809 1.845
# Sample mean method for comparison
set.seed(5)
r <- replicate(1000, mean(f(runif(n, a,b)))*(b-a))
hist(r)
summary(r)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1.788 1.798 1.800 1.800 1.803 1.813
Re your edit: I am assuming the x*2 + y^2, [-1,1] you are referring to a circle rather than a function f(z). So really to estimate area of unit circle/Pi by simulation.
f2 <- function(x) sqrt(1-x^2)
s <- seq(-1 , 1 ,by=0.001)
plot(s,f2(s))
# Get the max value of function within the range
c <- ceiling(max(f2(s)))
# [1] 1
n <- 1000000
a <- -1
b <- 1
set.seed(5)
x <- runif(n,a,b)
y <- c*runif(n,0,1)
R <- sum(y < f2(x))/n
(b-a)*c*R
#[1] 1.57063 # multiply it by 2 to get full area
pi/2
#[1] 1.570796