how can we generate random numbers using skew normal distribution in multivariate case how can we generate random numbers using skew normal distribution in multivariate case
If you are interested in theory, as @Nate Eldredge suggested, follow the url.

As written there also, you could use R in practice to generate random numbers using a skew normal distribution.

Load the library sn:

library(sn)


Or if not installed, first install it via install.packages('sn').

Then you can generate any random number with given parameters with the rsn function.

rsn(n=100, location=1.256269, scale=1.605681, shape=5)


Will generate 100 random numbers from the distribution with given parameters. If you would like to plot the histogram of your generated values, use higher sample size (e.g. 1.000 or 10.000), like:

hist(rsn(n=10000, location=1.256269, scale=1.605681, shape=5))
Edit: this answers a different question.

Supposing that the mean is $\mu$ and the variance is $\Sigma$, calculate the Cholesky decomposition $\Sigma =LL^*$, generate a vector $v$ of $\operatorname{rank} \Sigma$ independent Gaussians, and output $Lv + \mu$.