Generate random numbers in R with Skew-Normal Bivariate distribution

I made the script below, but I doubt if it is correct, I need to generate a bivariate vector with Skew-normal distribution of mean mu=(mu1,mu2) , variance and covariance matrix Sigma=matrix(c(sigma1^2,rho*sigma1*sigma2,rho*sigma1*sigma2,sigma2^2),2,2) and symmetry parameter lambda=(lambda1,lambda2) :

#Pacote utilizado
library("sn")
#Semente fixada
set.seed(0)
n=100 #Tamanho das amostras
#Valores iniciais dos parâmetros usados para gerar as variáveis. 
#Ou seja, Valor verdadeiro dos parâmetros
gamma0  <- 1
gamma1  <- 1
gamma2  <- 1.5
#gamma3 <- 0.7
beta0   <- 0.5
beta1   <- 1.5
#beta2  <- 1.1
sigma   <- 1
rho     <- 0.5
lambda1 <- 1
lambda2 <- 0
#Matriz de covariáveis para gerar mu1
X1 <- rep(1,n)
X2 <- rnorm(n,0,1)
X3 <- rnorm(n,0,1)
X <- cbind(X2)
#Matriz de covariáveis para gerar mu2
W <- cbind(X2,X3)
#Vetor de valores verdadeiros dos parâmetros
b <- rbind(beta0,beta1)
g <- rbind(gamma0,gamma1,gamma2)
#Vetor de médias 1
mu1 <- (as.numeric(model.matrix(~X) %*% b))
#Vetor de médias 2
mu2 <- (as.numeric(model.matrix(~W) %*% g))

Y <- matrix(NA,n,2)
  for(k in 1:n){
    Y[k,] <- rmsn(1,c(mu1[k],mu2[k]),matrix(c((sigma)^2,rho*sigma,rho*sigma,1),2,2),c(lambda1,lambda2))
  }

  Y1 <- cbind(Y[,1])
  Y2 <- cbind(Y[,2])

hist(Y1)
hist(Y2)

I am building a script for regression model analysis, so the averages are generated according to a linear model. Is my script correct?