Approximation of a multivariate function is an important theme in the field of numerical analysis and its applications, which continues to receive a constant attention. In this paper, we provide a parameterized family of self-referential (fractal) approximants for a given multivariate smooth function defined on an axis-aligned hyper-rectangle. Each element of this class preserves the smoothness of the original function and interpolates the original function at a prefixed gridded data set. As an application of this construction, we deduce a fractal methodology to approach a multivariate Hermite interpolation problem. This part of our paper extends the classical bivariate Hermite's interpolation formula by Ahlin (Math. Comp. 18:264-273, 1964) in a twofold sense: (i) records, in particular, a multivariate generalization of this bivariate interpolation theory; (ii) replaces the unicity of the Hermite interpolant with a parameterized family of self-referential Hermite interpolants which contains the multivariate analogue of Ahlin's interpolant as a particular case. Some related aspects including the approximation by multivariate self-referential functions preserving Popoviciu convexity are given, too.