A new method is presented for simulating discrete Coulomb interaction effects in electron and ion beam columns. The method is applicable to the latest high-throughput electron and ion beam projection systems currently under development as possible candidates for next generation lithography. Monte Carlo simulation is used to trace bunches of several thousand particles simultaneously down the column. Compared with previous simulations, two key improvements have been made: (1) Instead of using thin lens approximations, the real fields of the magnetic and electrostatic lenses are used, which enables the lens properties and aberrations to be accurately included in the computation; (2) The N-body inter-particle Coulomb forces are computed with a hierarchical tree-code algorithm, which is orders of magnitude faster that the direct pair-wise force evaluation method used previously. The real lens fields are accurately computed with second-order finite element method ('SOFEM'), and the lens fields are then fitted with Fourier-Bessel series that accurately represent the real lens fields while simultaneously being exact solutions of Laplace's equation. The discrete inter-particle Coulomb fields, evaluated with the new tree-code algorithm, are then added to the external fields of the lenses. Trajectories through these combined fields are then computed by direct ray-tracing, using a fifth-order Runge-Kutta formula, to an internal self- consistency of better than 1 picometer. This method enables, for the first time, the combined effects of the discrete Coulomb interactions and the lens aberrations to be predicted accurately in a completely unified and self-consistent way. The method is illustrated for the case of a magnet doublet projection system. The method can equally well be applied to multi-beam columns, electrostatic lenses, cathode lenses and electron mirrors.