Within a cosmological context, we study the behaviour of collisionless particles in the weak field approximation to General Relativity, allowing for large gradients of the fields and relativistic velocities for the particles. We consider a spherically symmetric setup such that high resolution simulations are possible with minimal computational resources. We test our formalism by comparing it to two exact solutions: the Schwarzschild solution and the Lemaître-Tolman-Bondi model. In order to make the comparison we consider redshifts and lensing angles of photons passing through the simulation. These are both observable quantities and hence are gauge independent. We demonstrate that our scheme is more accurate than a Newtonian scheme, correctly reproducing the leading-order post-Newtonian correction. In addition, our setup is able to handle shell-crossings, which is not possible within a fluid model. Furthermore, by introducing angular momentum, we find configurations corresponding to bound objects which may prove useful for numerical studies of the effects of modified gravity, dynamical dark energy models or even compact bound objects within General Relativity.