Orbital propagation becomes a High-Performance Computing (HPC) problem when many trajectories need to be simulated, e.g., for space debris modelling or planetary protection analysis. In this work, Alessandro Masat, Camilla Colombo, and Arnaud Boutonnet (ESA-ESOC) investigate the augmentation of the Picard-Chebyshev numerical scheme to maximise the computational performance of the propagation of large sets of initial conditions. The simulation method is studied both on CPU-based and GPU-based computational architectures, with the latter that achieves server-like performance using a single low-end graphics card only.