Abstract A novel computational scheme that models the scattering of light at interfaces in a deep‐nanometric (deep‐nm) plasmonic system is presented. The nonclassical effects associated with the collective motion of free electrons in metals of the system are investigated by a nonlocal hydrodynamic (HD) model. Due to the HD model, three types of interfaces, that is, the dielectric–dielectric interface, the dielectric–metal interface, and the metal–metal interface, are distinguished. The scattering of light at these interfaces is described by Boundary Integral Equations (BIEs), which are constructed by using the surface equivalence principle, the concept of the Green' function and the interface conditions. The BIEs are numerically solved by the boundary element method (BEM) for 3D nanoparticles of arbitrary shapes. To validate, the computed results are first quantitatively contrasted with analytical solutions from the Mie theory for spherical core–shell structures, with good agreements being demonstrated in both the near field and the far field. Then, a qualitative study is performed for dimers consisting of spheres with various gap sizes for both the classical local response model and the nonlocal HD model. The observed trends agree well with the results reported in previous works.