A Monte Carlo program is presented which simulates the response of SiPMs in the nonlinear regime, where the number of Geiger discharges from photons and/or from dark counts in the time interval given by the pulse shape of a single Geiger discharge, approaches or exceeds the number of SiPM pixels. The model includes the effects of after-pulses, of prompt and delayed cross-talk, and of the voltage drop over the load resistance of the readout electronics. The results of the simulation program are compared to published results from SiPMs with different number of pixels for different intensities and time distributions of photons, dark-count rates, SiPM pulse shapes, and probabilities of cross-talk and after-pulsing.