A procedure is derived for simulating correlated non-normal systems of linear statistical equations. The method is based on fifth-order polynomial transformations to generate multivariate non-normal distributions. The procedure allows for the simultaneous control of the correlated non-normal (a) stochastic disturbance distributions, (b) independent variables, and (c) dependent and independent variables for each equation throughout a system. A numerical example is provided to demonstrate the procedure. The results of a Monte Carlo simulation are provided to confirm that the proposed method generates the specified standardized cumulants and correlations.