A simple model to simulate cycle-by-cycle variation that is suitable for use in Monte-Carlo approaches has been developed and validated with a wide range of experimental data. The model is intended to be diagnostic rather than predictive in nature, with a goal of providing realistic in-cylinder pressures. The individual-cycle cumulative rate of heat release was curve fit with a four-parameter Wiebe function. It was found that the distribution of the Wiebe b-parameter was quite small, so its value was obtained from the ensemble-averaged condition. The remaining three Wiebe function parameters, θig, θcomb and m were found to be distributed over a moderate range, and were linearly correlated to each other. Using the cumulative density function of θig, and the linear fit of θcomb and m to θig, with a random component added, a Monte-Carlo scheme was developed. The resulting randomly chosen cycles were found to adequately reproduce the observed cyclic variation at the same condition, and can be used with higher-order models to assess the impact of cyclic variations on other cylinder-pressure-driven phenomena.