Dependencies among input variables in Monte-Carlo simulation studies constitute a well-known problem in computational statistics. Dependencies between random variables are usually described either by Pearson's coefficients of linear correlation or by Spearman's rank-order correlations. The predominant use of the latter is a consequence of the widely held belief that it is too difficult to maintain a specified Pearson correlation matrix for other than normally distributed random variables. On the contrary, there exists a procedure by Li and Hammond (IEEE Trans. Systems Man Cybernet 5 (1975) 557-561) that exactly achieves this. Their solution is based on a translation method in which the standard multi-variate normal distribution is transformed into the arbitrary multi-variate distribution. It requires the numerical computation of a key integral, involving transformations to the marginal distributions, that becomes singular at the extreme correlations. However, if following Lancaster (Biometrika 44 (1957) 289-292) the transformations are expanded in Chebyshev-Hermite polynomials, a simple explicit expression of the integral expansion coefficients results. This paper stabilizes the actual computations using the recurrence relation for Chebyshev-Hermite polynomials. The new solution is fast and can be made arbitrarily accurate. An implementation for piecewise linear marginal distributions that handles continuous, discrete and mixed distributions covers most practical purposes. A number of examples show the potential of the new method. (~) 1998 P.A.G. van der Geest. Published by Elsevier Science B.V.