A singularity method is presented for solving the linearized Poisson-Boltzmann equation in order to determine electrostatic interactions in multiparticle systems, Results are presented for (1) two and three unequal spheres interacting in bulk solution, (2) a single sphere interacting with a planar interface, (3) a sphere interacting with a planar interface on which three other spheres are adsorbed, and (4) a sphere interacting with a cylindrical fiber or rod-like particle, Boundary conditions of constant surface potential, constant surface charge, and the "linearized regulation model" are used, The results obtained are shown to be in excellent agreement with existing approximations where appropriate and with existing exact solutions where available, Calculations involving more than two spheres are compared with results derived by using pairwise additive interactions, and limitations of the pairwise additivity approximation are discussed, The method is developed in a way that is not specific to any particular locations of the particles or boundaries, making it possible to incorporate it into Stokesian dynamics or Monte Carlo simulations of many interacting spheres.