library(MASS)

# set correlation coefficient and make matrix
r <- 0.9
(corr_matrix <- matrix(c(1, r, r, 1), nrow = 2))
##      [,1] [,2]
## [1,]  1.0  0.9
## [2,]  0.9  1.0
# generate multivariate normal distribution
tmp <- mvrnorm(
  n = 5000,
  mu = c(0,0),
  Sigma = corr_matrix
)

# get probability that random variable with N(0,1) is below the value given
dat <- pnorm(tmp)

# check that the variables are uniformly distributed between 0 and 1
par(mfrow = c(1,2))
hist(dat[,1])
hist(dat[,2])

# check that the variables are correlated
par(mfrow = c(1,1))
plot(dat[,1], dat[,2], pch = 20, cex = 0.05)