Shepard’s method is a fast algorithm that has been classically used to interpolate scattered data in several dimensions. This is an important and well-known technique in numerical analysis founded in the main idea that data that is far away from the approximation point should contribute less to the resulting approximation. Approximating piecewise smooth functions in \(\mathbb {R}^{\varvec{n}}\) near discontinuities along a hypersurface in \(\mathbb {R}^{\varvec{n-1}}\) is challenging for the Shepard’s method or any other linear technique for sparse data due to the inherent difficulty in accurately capturing sharp transitions. This article is devoted to constructing a non-linear Shepard’s method using the basic ideas that arise from the weighted essentially non-oscillatory interpolation method (WENO). The proposed method aims to enhance the accuracy and reduce the smearing of the traditional Shepard’s method by incorporating WENO’s adaptive and non-linear weighting mechanism. To address this challenge, we non-linearly modify the weight function in a general Shepard’s method, considering any weight function, rather than relying solely on the inverse of the distance squared. This approach effectively reduces the smearing of discontinuities providing a sharper approximation. The numerical experiments presented demonstrate the superior performance of the new method close to the discontinuities and confirm the theoretical results exposed in this manuscript.