We introduce an efficient and accurate boundary element method for computing the electrostatic potential in closely-packed dielectric spheres and Janus particles. Close physical proximity occurs in modeling molecular ion interactions or colloidal materials. The electrostatic potential, which is described by the Poisson-Boltzmann equation, becomes highly singular under close interactions, resulting in difficulties in computational results. We select a boundary element approach due to the natural reduction of the dimension of the problem from volume to surface elements. The singular behavior caused by the close proximity of the system due to the source charge and part of the reaction field is removed from the potential via a subtraction de-singularization technique within a hyper-singular, high order second kind integral equation formulation. The resulting system of equations has a number of right-hand-side integrals that do not contain basis elements and contain the bulk of the singular behavior. These auxiliary integrals require treatment in order to best capture the singular behavior while keeping cost at a minimum. So long as these integrals are calculated accurately, the system can be readily, easily solved. Regularization techniques for the Hadamard finite part integral that appears in this method are then presented, where mathematical identities and adaptive meshes offer a means to compute the singular integral with the required level of accuracy at a much reduced computational cost. Adaptive quadrature is presented and used for other less challenging, but still singular integrals. We demonstrate that numerical results of the singular potential for one and two closely-packed spheres have validated the effectiveness and accuracy of the proposed method.