Multislice simulations of 4D scanning transmission electron microscopy (4D STEM) data are computationally demanding due to the large number of STEM probe positions that must be calculated. For accurate analysis, inelastic scattering from phonons and plasmons must also be included. However, current frozen phonon and Monte Carlo plasmon techniques require a separate calculation for each different phonon/plasmon configuration, and are therefore not suitable for scaling up to 4D STEM. Here a phase scrambling algorithm (PSA) is proposed, which treats all phonon/plasmon configurations simultaneously. A random phase is introduced to maintain incoherence between the different inelastic scattering events; this is the phase scrambling part of the algorithm. While for most applications, a few tens of frozen phonon iterations are sufficient for convergence, in the case of plasmon scattering as many as tens of thousands of iterations may be required. A PSA is statistically more representative of inelastic scattering, and achieves significant savings in computation time for plasmons. The increase in speed is a pre-requisite for 4D STEM inelastic scattering simulations.