library(MASS)

set.seed(142857)
mu_vec <- c(0, 0, 0)  # The means
cov_mat <- diag(1, 3) # The covariance matrix
n <- 100

x <- mvrnorm(n = n, mu = mu_vec, Sigma = cov_mat, empirical = TRUE)

cov(x)
              [,1]          [,2]          [,3]
[1,]  1.000000e+00 -2.250771e-16 -2.604286e-16
[2,] -2.250771e-16  1.000000e+00  6.497218e-17
[3,] -2.604286e-16  6.497218e-17  1.000000e+00