Hi,
The Statistics Toolbox has a function mvnrnd to generate jointly Gaussian random variables with specified means and covariance matrix:
>>N = 10; % desired number of samples of each variable
>>mu = [10; 10]; % vector of means
>>cov = [3 1; 1 3]; % covariance matrix
>>samples = mvnrnd(mu, cov, N);
If you want to do it manually, you can generate independent standard Gaussian RV's (with randn) and apply an affine transformation that will give the desired mean vector and covariance matrix. For that you need to compute the Cholesy decomposition of the latter. Note that the output of Matlab's chol function (see its documentation) needs to be (conjugate-)transposed to conform to the procedure described in the link.
The code would be:
>>N = 10; % desired number of samples of each variable
>>mu = [10; 10]; % vector of means
>>cov = [3 1; 1 3]; % covariance matrix
>>independent_samples = randn(N, size(cov,1));
>>A = chol(cov)';
>>samples = bsxfun(@plus, mu, A*independent_samples.').';
