We present a computational scheme based on classical molecular dynamics to study chaotic billiards in static external magnetic fields. The method allows us to treat arbitrary geometries and several interacting particles. We test the scheme for rectangular single-particle billiards in magnetic fields and find a sequence of regularity islands at integer aspect ratios. In the case of two Coulomb-interacting particles the dynamics is dominated by chaotic behavior. However, signatures of quasiperiodicity can be identified at weak interactions, as well as regular trajectories at strong magnetic fields. Our scheme provides a promising tool to monitor the classical limit of many-electron semiconductor nanostructures and transport systems up to high magnetic fields.