I have a simple code that performs a particle simulation here: 2021_FortranCon/particle_simulation.f90 at main · m3g/2021_FortranCon · GitHub
With periodic boundary conditions. It ends up producing this kind of animation, when visualized with vmd: Simple 2D particle simulation - YouTube
(there is a script in the same repository to visualize the trajectory generated with the Visual Molecular Dynamics (VMD) package).
This kind of animation of particles bouncing usually traps the attention of students. It is easy to adapt the code to perform simulations of planetary motion in a gravitational field (by removing the periodic boundary conditions).