The simulation of particle-laden flows is a crucial task in fluid dynamics, requiring high computational cost owing to the complex interactions between numerous particles. Typically, the flow velocity is described with the equations proposed by Stokes. While there is an analytical solution for the Stokes flows around a single spherical particle, the Stokes flows around many particles are still unsolved. In this paper, we study Genetic Programming (GP) for symbolic regressions to explore the potentials of multi-objective GP in recovering analytical expressions for two and, in the future, N particles. We propose a new GP approach containing building blocks to scale up the problem and provide a new benchmark with 22 cases for this application. To identify the strengths and limitations of GP, we generate fully resolved training data from simulations. We compare the results of our algorithm to the superimposition method and a multi-layer perceptron as two baseline methods. The results show that GP can find comparable and sometimes better solutions with smaller failure rates than the two baseline methods. In addition, the produced solutions by GP are explainable and certain function patterns inline with physical laws can be identified across the benchmark problems.