Extending of given sequances

Lets say, we have two given sequences X and Y. How to extend these sequences that them mean, dispersion, correlation coefficient and other parameters equivalent. If I need only need the mean and dispersion equivalent, I can generate X2 = N(0,1) * sqrt(DX) + EX, and Y2 = N(0,1) * sqrt(DY) + EY, where N - normal distribution, D - dispersion, E - mean. But how should i generate X2 and Y2, that their correlation coefficient do not vary?
I tried to figure out an algorithm. Lets say, we have figured out X and Y dependence equation X = c1 * Y + c2. So, we can generate Y2 = N(0,1) * sqrt(DY) + EY and X2 by given dependence, but we have to modify it - we need some variation about correlation. We add another value R = (1-corr) + Rand() * (1+corr-1+corr) = (1-corr) + 2*corr * Rand(), where corr - correlation coefficient, when R mean is 1. So X2 = (c1 * Y2 + c2) * R, which mean should be equal EX (but what about DX?). So I extended X and Y with many generated new values and checked the results. The means is ok, but X dispersion was a little different from given and new correlation coefficient is similar to given, but a little different too.
So, what do I do wrong? I formulated that algorithm by myself, so may there is another real algorithm? Thanks for help.

To generate random samples from the bivariate normal distribution, do the following. Generate \(\displaystyle Z_1, Z_2\) as standard normal variates. Then let \(\displaystyle X = \mu_1 + \sigma_1 Z_1, Y = \mu_2 + \sigma_2 ( Z_1 \rho + Z_2 \sqrt{1-\rho^2} })\).

Thanks a lot