The orbital propagation of large sets of initial conditions under high accuracy requirements is currently a bottleneck in the development of space missions, e.g. for planetary protection compliance analyses. The proposed approach can include any force source in the dynamical model through efficient Picard-Chebyshev (PC) numerical simulations. A two-level augmentation of the integration scheme is proposed, to run an arbitrary number of simulations within the same algorithm call, fully exploiting high performance and GPU (Graphics Processing Units) computing facilities. The performances obtained with implementation in C and NVIDIA CUDA programming languages are shown, on a test case taken from the optimization of a Solar Orbiter-like first resonant phase with Venus.