This article is the second in a series of two papers concerning the mathematical study of a boundary integral equation of the second kind that describes the interaction of N dielectric spherical particles undergoing mutual polarisation. The first article presented the numerical analysis of the Galerkin method used to solve this boundary integral equation and derived N-independent convergence rates for the induced surface charges and total electrostatic energy. The current article will focus on computational aspects of the algorithm. We provide a convergence analysis of the iterative method used to solve the underlying linear system and show that the number of liner solver iterations required to obtain a solution is independent of N. Additionally, we present two linear scaling solution strategies for the computation of the approximate induced surface charges. Finally, we consider a series of numerical experiments designed to validate our theoretical results and explore the dependence of the numerical errors and computational cost of solving the underlying linear system on different system parameters.