The focus on the covariance matrix seems to imply that the data are drawn from a Gaussian distribution. So a starting point would be the functions of the GNU scientific library, https://www.gnu.org/software/gsl/manual/html_node/The-Gaussian-Distribution.html#The-Gaussian-Distribution.

If you're using a statistics package, it probably has a library function to do this. R, Matlab, and Mathematica all do, assuming you want to generate a multivariate Gaussian with the required covariance. Otherwise, generate vectors from an isotropic Gaussian with unit variance and multiply them by one of the matrices of a Cholesky factorisation of the covariance matrix.

Assuming you need two sets of data that has prescribed correlation. If the correlation matrix R that is positive definite is prepared. and the two normally distributed parameters are sets of monte carlo realizations Wi are prepared than one of the possibility is via application of Choleski decomposition Q=chol(R) and set of correlated normal distributions Yi are Yi = Wi*R.

The procedure is described through equations (5)-(8) in linked paper. I think that it will be clear even though it is in Czech. It also gives procedure how to trasform to different than normally distributed parameters.

Petr

Article Simulace korelovaných neparametrických rozdělení v rámci metody SBRA