A particle-source-in-cell (PSIC) method is developed that is based on a high-order discontinuous Galerkin method for the simulation of two-dimensional particle-laden flows with a mutual (two-way) coupling between an inertial particle phase and the carrier flow phase on parallel computer architectures. Starting from a parallel, high-order discontinuous Galerkin (DG) code for the inviscid Euler equations with a one-way coupled parallel particle tracer, a smooth polynomial shape function is implemented that distributes each particle’s forcing influence on the Eulerian grid with high-order accuracy. A parallel computing capability is developed that is consistent with the existing particle tracking algorithm through MPI communication protocols. In a pre-computation procedure, the non-overlapping regions that identify the nearest node in the mesh are established first. Then the 4-by-4 element span connected to this node is determined in order to ensure the sufficient support for the compactly-supported particle distribution function. This information is stored on each processor and is used to distribute the particle’s influence efficiently on parallel computers. After verification of the shock capturing capability of the Euler code against the classical Sod one-dimensional shock tube test case, the PSIC solver is tested for two one-dimensional flow cases, including (1) a particle cloud accelerated in a uniform subsonic flow and (2) a particle cloud accelerated in the supersonic flow behind a moving shock wave. These computations are performed using a rectangular domain with a structured grid of orthogonal elements. Periodic boundary conditions are specified in the cross-stream direction so that the flow is nominally one-dimensional. The time-dependent pressure profiles and number density profiles along the streamwise coordinate for both test cases are in excellent comparison against published data of PSIC computations based on a Weighted Essentially Non-Oscillatory (WENO) finite difference method. The results show no discernible variation in the cross-stream coordinate which validates the high-order DG-based PSIC solver.