Several methods have been posted; here's another.
If you can sample from the $2$-dimensional normal distribution whose variance is $\begin{bmatrix} 1 & \rho \\ \rho & 1\end{bmatrix}$, you can then rescale the $x$-coordinate and the $y$-coordinate to get the standard deviations you want, and that won't alter the correlation $\rho$. If $A$ is a square root of the given correlation matrix and $\begin{bmatrix} Z_1 \\ Z_1 \end{bmatrix}$ has a $2$-dimensional normal distribution whose variance is the $2\times 2$ identity matrix, then $A\begin{bmatrix} Z_1 \\ Z_1 \end{bmatrix}$ has a $2$-dimensional normal distribution whose variance is $\begin{bmatrix} 1 & \rho \\ \rho & 1\end{bmatrix}$.
So the remaining problem is finding a square root of $\begin{bmatrix} 1 & \rho \\ \rho & 1\end{bmatrix}$.
Let $P=\begin{bmatrix} 1/2 & 1/2 \\ 1/2 & 1/2\end{bmatrix}$ and $Q=I-P$. Then $P^2=P$, $Q^2=Q$, and $PQ=QP=0$.
Then we have $\begin{bmatrix} 1 & \rho \\ \rho & 1\end{bmatrix} = (1+\rho)P + (1-\rho)Q.$ Since $P^2=P$, $Q^2=Q$, and $PQ=QP=0$, we can find a square root just by taking square roots of the coefficients: $ \sqrt{1+\rho}\,P + \sqrt{1-\rho}\,Q. $