This paper presents a highly efficient and stable algorithm for the numerical simulation of collision free plasma. The algorithm has been successfully used to numerically model non linear electron cyclotron resonance in VLF band radio waves in space, and has produced good simulations of radio emissions such as ‘dawn chorus’ and ‘triggered VLF emissions’. The algorithm fills the phase box with simulation particles which represent phase space trajectories. Particle trajectories are followed forwards continuously in time without restarts, and these particles are followed until either the simulation ends or the particles exit from the phase box. Liouville’s theorem states that distribution function F is conserved along these trajectories, and so F is known at the phase space points represented by the particle positions. No phase space volume is associated with these particles. At each time step F is truly interpolated, and not distributed from particles to the regular phase space grid, solely for the purpose of calculating plasma charge and current fields. For this operation a cheap low order interpolation is useable. The algorithm is simple, robust and highly efficient. There is NO diffusion in phase space whatsoever, and fine structure is merely undersampled and does not cause any instability problems. Regions of phase space where F=0 are cost free as no particles need be provided. The algorithm readily accomodates situations where phase fluid flows into or out of the phase box. The paper presents the latest simulations of VLF chorus and VLF triggered emissions from the VHS/VLF code, which are in excellent agreement with observations.