In this dissertation we develop a high-order regularization technique with optimal scaling to approximate singular sources expressed as one single Dirac-delta or weighted summation of Dirac-deltas. We consider the numerical solution of hyperbolic conservation laws arising in the simulation of shocked particle-laden flows with the particle-source-in-cell (PSIC) method. In one dimension, the regularization technique is based on a class of compactly-supported piecewise polynomials plus an optimal scaling of the width of the support. We show that the regularization approximates the singular source with the desired order of accuracy and smoothness (number of continuous derivatives). We establish a theoretical criterion to choose an optimal scaling in the regularization that leads to formal order of convergence of the numerical scheme, away from the singularities, in the solution of the hyperbolic conservation laws. We validate our criterion by solving a linear and a nonlinear (Burgers) scalar hyperbolic conservation law with a singular source, as well as the nonlinear Euler equations with singular sources, a system of hyperbolic conservation laws governing compressible fluid dynamics with shocks and particles. A Chebyshev collocation method (spectral) discretizes the spatial derivatives in the scalar equation tests. A multidomain high-order resolution hybrid spectral-WENO method discretizes the Euler equations. In the future, we aim to enhance the accuracy of the PSIC method by using the proposed regularization technique with optimal scaling.