Given a $M \times M$ desired covariance, $R$, and a desired number of sample vectors, $N$ calculate a $N \times M$ Gaussian random vector, $X$.
Not really sure what to do here. You can calculate the joint pdf given a mean, $\mu$, and covariance. So for $2 \times 2$ covariance matrix is defined as: $$ \text{Cov}[X] = \begin{bmatrix}\text{Cov}(X_1) & \text{Cov}(X_1,X_2) \\ \text{Cov}(X_2,X_1) & \text{Cov}(X_2)\end{bmatrix} $$ But not sure how to get the mean from that.
To generate one vector $u\in\mathcal{R}^M$, first of all generate any vector $v$ from $\mathcal{N}(0,I)$ (or M independent normally distributed variables with mean $0$, varaince $1$).
We now need to get a matrix $L$ such that $LL^T=R$, easiest way is a cholesky decomposition. Now $u=Lv\sim\mathcal{N}(0,R)$