One suggestion is to work with copulas. In a nutshell, a copula allows you to separate out the dependency structure of a distribution function. Say, $F_1,F_2,\ldots,F_n$ are the 1D marginals of a distribution $F$ then the copula $C$ is the function defined as
$$C(u_1,u_2,\ldots,u_n)=F(F^{-1}_1(u_1),F^{-1}_1(u_2),\ldots,F^{-1}_n(u_n))$$
This makes $C$ a function from $[0,1]^n$ to $[0,1]$. For instance, if you take the bivariate normal distribution, by doing the computation above, you'll find the Gaussian copula
$$C^{\text{Gauss}}_{\rho}=\int_{-\infty}^{\phi^{-1}(u_1)}\int_{-\infty}^{\phi^{-1}(u_2)}\frac{1}{2\pi\sqrt{1-\rho^2}}\exp\left(-\frac{s_1^2-2\rho s_1s_2+s_2^2}{2(1-\rho^2)}\right)ds_1ds_2$$
I used the package copula in R to illustrate. If you just take the copula as such, it is as if you constructed a probability distribution with the dependency structure of a bivariate normal, but with uniform marginals. So I generated 1000 random vectors from a Gaussian copula with $\rho=0.5$. Here's the code
library(copula);
norm.cop <- normalCopula(0.5);
u <- rCopula(1000, norm.cop);
plot(u,col='blue',main='Random variables, uniform marginals, gaussian copula, > rho=0.5',xlab='X1',ylab='X2')
cor(u);
and the result 
I also computed the sample correlation which is $0.5060224$.
I also computed a plot to show you the marginals are indeed uniform
dom<-(1:length(u[,1]))/length(u[,1]);
par(mfrow=c(1,2));
plot(dom,sort(u[,1]),col='blue',main='marginal X1');
abline(0,1,col='red');
plot(dom,sort(u[,2]),col='blue',main='marginal X2');
abline(0,1,col='red');
This is all very nice, but there are a number of pitfalls that have to be discussed:
- Copula's for discrete distributions are a real can of worms.
- If we can use a multivariate Gaussian distribution to get a dependency structure, why not use a multivariate student t? Or a multivariate Pareto? Or other dependencies which are much more exotic, but all could in principle also lead to a $0.5$ correlation if you set the parameters right.
- Given marginals and a correlation, it is not always the case that you can construct a copula and hence a multivariate distribution with the desired properties. A nice example is given in Embrechts (2009), "Copulas: A Personal View", The Journal of Risk and Insurance, Vol. 76, No. 3, 639-650. The example shows that if you want the marginals to be lognormal distributed $LN(0,1)$ and $LN(0,6)$ respectively, your correlation is restricted to the range $[-0.00025,0.01372]$. The heavy tails of the lognormals essentially constrain you to that range. This can be proven from the Fréchet-Hoeffding bounds. More details are in the article.
More can be said and I think the article I quoted in my last item is a nice starting point.
Best Answer
It is a relatively simple task to generate samples of random variables with given marginal distributions that are correlated. The difficulty lies in controlling the exact degree of correlation, if that is desired, unless the marginal distributions are normal.
The Cholesky approach mentioned works well for constructing random variables with a multivariate normal distribution and a specified correlation matrix given a set of independent random variables with normal marginal distributions. For example suppose independent random variables $Z_1$ and $Z_2$ both have standard normal marginal distributions, i.e. $Z_1, Z_2 \sim N(0,1)$, then take
$$X = Z_1, \,\,\, Y = \rho Z_1 + \sqrt{1 - \rho^2}Z_2.$$
Such a transformation preserves the marginal standard normal distributions, i.e. $X, Y \sim N(0,1)$ and imposes the desired correlation
$$E(XY) = \rho E(Z_1^2) + \sqrt{1- \rho^2}E(Z_1Z_2) = \rho.$$
An approximate approach for non-normal marginal distributions, $F$ and $G$, would be to first draw independent samples from a standard normal distribution, $Z_1, Z_2 \sim N(0,1)$. Next impose a correlation $\rho$ using the transformation
$$V_1 = Z_1, \,\,\, V_2 = \rho Z_1 + \sqrt{1-\rho^2}Z_2.$$
Note that $V_1$ and $V_2$ have a joint normal distribution. If $\Phi$ is the standard normal cummulative distribution function then $\Phi(V_1)$ and $\Phi(V_2)$ have uniform $U(0,1)$ distributions, since, for example,
$$P(\Phi(V_1) \leqslant v) = P(V_1 \leqslant \Phi^{-1}(v)) = \Phi[\Phi^{-1}(v)] = v. $$
Finally perform the following transformation using inverse marginal distribution functions $F^{-1}$ and $G^{-1}$ and the standard normal cumulative distribution function $\Phi$,
$$X = F^{-1}[\Phi(V_1)], \,\,\, Y = G^{-1}[\Phi(V_2)].$$
Now $X$ and $Y$ have the desired marginal distributions since, for example,
$$P(X \leqslant x) = P(F^{-1}[\Phi(V_1)] \leqslant x) = P(\Phi(V_1) \leqslant F(x)) = F(x).$$
In general due to non-linearity, $corr(X,Y) \neq \rho$, but it may not be far off and you can iterate on the choice of $\rho$ in the first step until you get close to the desired correlation.
A more comprehensive treatment of imposing a dependence structure on random variables with given marginals can be found in the theory of copulas.