Simulating multivariate nonnormal data with specified correlation matrices is difficult. One especially popular method is Vale and Maurelli's (1983) extension of Fleishman's (1978) polynomial transformation technique to multivariate applications. This requires the specification of distributional moments and the calculation of an intermediate correlation matrix such that when data are transformed, the target correlation matrix is reproduced. We present an alternative technique that involves sampling data from specified population distributions and identifying the intermediate correlation matrix through an iterative, trial-and-error process. We provide R program code to implement this technique and show that it can generate data under a wide range of conditions (e.g., with empirical samples, with discrete rather than continuous data, when distributional moments are undefined or outside the boundary conditions of other techniques). This approach could be useful in many contexts, especially Monte Carlo studies of multivariate statistics. (Contains 5 tables, 3 figures and 2 footnotes.)