In this work, we develop and analyze a formalism for solving boundary value problems in arbitrarily-shaped domains using the MADNESS (multiresolution adaptive numerical environment for scientific simulation) package for adaptive computation with multiresolution algorithms. We begin by implementing a previously-reported diffuse domain approximation for embedding the domain of interest into a larger domain (Li et al., 2009 [1]). Numerical and analytical tests both demonstrate that this approximation yields non-physical solutions with zero first and second derivatives at the boundary. This excessive smoothness leads to large numerical cancellation and confounds the dynamically-adaptive, multiresolution algorithms inside MADNESS. We thus generalize the diffuse domain approximation, producing a formalism that demonstrates first-order convergence in both near- and far-field errors. We finally apply our formalism to an electrostatics problem from nanoscience with characteristic length scales ranging from 0.0001 to 300 nm.