Verifying if a function is 2-increasing

Question

In the book An Introduction to Copulas (2006), by Roger B. Nelsen, on page 8, the author defines a function $H$ as 2-increasing if

\begin{equation} V_H(B)=H(x_2,y_2)-H(x_2,y_1)- H(x_1,y_2)+H(x_1,y_1)\geq0 \end{equation} for all rectangles $B$ whose vertices lie in Dom$H$.

Where $V_H(B)$ is called de $H$-volume 0f $B$. The rectangle $B$ in $\overline{\mathbb{R}}\times\overline{\mathbb{R}}$ is the Cartesian product of two closed intervals $[x_1,x_2]\times[y_1,y_2]$ whose vertices are the points $(x_1,y_1),(x_1,y_2),(x_2,y_1),(x_2,y_2)$. $\overline{\mathbb{R}}$ is the extended line of real numbers $[-\infty;\infty]$.

Later in the book, on pages 17-18, it is said that the function

H(x,y)= \begin{cases} \frac{(x+1)(e^y-1)}{x + 2e^y -1}, & (x,y)\in [-1,1]\times[0,\infty] \\[2ex] 1-exp(-y), & (x,y)\in (1,\infty]\times[0,\infty]\\[2ex] 0, & elsewhere. \end{cases}

is 2-increasing.

I wanted to verify this statement and I decided to simulate $V_H(B)$ in R for the part $H(x,y)=\frac{(x+1)(e^y-1)}{x + 2e^y -1}$ to check if for every computation $V_H(B)\geq 0$.

After one million simulations, I have many NaN and negative results. Does it mean that the function is not 2-increasing ? If the function is 2-increasing, where am I mistaken ?

Thank you in advance.

My R code :

H<- function(x, y){

  ((x+1)*((exp(y))-1))/(x + 2*exp(y) -1) #Function H(x,y)
}

Hxy <- c()

for(i in (1:1000000)){

  x <-sort(runif(2,-1,1)) # generate randomly and sort
                          # x1 and x2 at each iteration

  y <-sort(runif(2,0,1000)) # generate randomly and y1 and y2 
                            # at each iteration
  
  Hxy[i] <- H(x[2],y[2])- H(x[2],y[1]) - H(x[1],y[2]) + H(x[1],y[1]) 
            # Compute H-volume of B
  
  i = i + 1
}