The proliferation of experimental platforms at the deep-nanometric scale in the last two decades has allowed to put into test theoretical results stemming from classical electromagnetics and pertaining to the vibrant field of nanoplasmonics. Ever-increasing deviations with the scaling-down of dimensions between experiment and classical results caused a renewed interest in semiclassical models, prominently among them the Hydrodynamic Drude Model (HDM). HDM and similar approaches mimic reasonably accurately intrinsically microscopic phenomena (e.g., nonlocal electron dynamics for the HDM) and in a computationally appealing format in comparison with full ab initio theories. Yet, their full computational potential can only be unlocked through numerical recipes that are particularly tailored to the new multiphysics and multiscale problems and which acknowledge the new physics that emanates from semiclassical approaches. We present here developments on a potential-based Boundary Integral Equation (BIE) method. This surface method bypasses issues due to the multiscale nature of the problem and by introducing a potential-based formalism, it accounts for the new degrees of freedom, which arise due to nonlocality in a lenient manner. We validate our recipe for the case of a simple sphere and showcase its application for a common theoretical platform for nanoplasmonics, the spherical dimer.