Phoretic particles exploit local self-generated physico-chemical gradients to achieve self-propulsion at the micron scale. The collective dynamics of a large number of such particles is currently the focus of intense research efforts, both from a physical perspective to understand the precise mechanisms of the interactions and their respective roles, as well as from an experimental point of view to explain the observations of complex dynamics as well as formation of coherent large-scale structures. However, an exact modelling of such multi-particle problems is difficult and most efforts so far rely on the superposition of far-field approximations for each particle's signature, which are only valid asymptotically in the dilute suspension limit. A systematic and unified analytical framework based on the classical Method of Reflections (MoR) is developed here for both Laplace and Stokes' problems to obtain the higher-order interactions and the resulting velocities of multiple phoretic particles, up to any order of accuracy in the radius-to-distance ratio ε of the particles. Beyond simple pairwise chemical or hydrodynamic interactions, this model allows us to account for the generic chemo-hydrodynamic couplings as well as N -particle interactions (N ≥ 3). The ε 5 -accurate interaction velocities are then explicitly obtained and the resulting implementation of this MoR model is discussed and validated quantitatively against exact solutions of a few canonical problems.