Management and exploitation in mines require a continuous and relatively smooth surface of the mineral grades. While assessing the various mineral elements, the scattered exploratory cavities are irregularly excavated. Producing a continuous surface from measured data requires interpolation methods. Several factors, including the characteristics of the data, affect the efficiency of the interpolation methods. For this reason, the efficiency of different methods in various cases is inconsistence, and choosing the appropriate interpolation method is also challenging. Interpolation methods can be categorized into two groups of mesh-based and meshless methods. Despite the efficiency and capabilities of meshless methods, they have a fundamental shortcoming due to the fixed size of the support domain. On the one hand, the distribution of exploratory cavities in mines is usually irregular, and in some areas, it is very dense, and in others, it is very sparse. On the other hand, the grade values of minerals at the surface of the region can be very variable with high changes. Conventional interpolation methods do not have sufficient efficiency and flexibility in confronting these two aforementioned issues. In this study, a precise, reliable, and flexible method is developed for interpolation of minerals through integrating the moving least squares and recursive least squares methods. In the proposed method for crack detection, the residuals statistical test of least squares computations is used.  In this method, for the central point, a continuity threshold (non-continuity) is determined based on the standard deviation of field values, so that points with crack are revealed and removed from the calculation of the value of the central point. Moreover, the size of the support domain is determined dynamically based on the recursive property of the method. In this method, an individual radius for the support domain is assigned to each central point according to the values and distributions of the surrounding field points. The dynamic size of the support domain allows a precise and reliable estimation of polynomial coefficients and the values of the central points. The efficiency of the proposed method is evaluated by applying it to simulated data as well as comparing it with the results of conventional interpolation methods on real mineral data. The results of the simulation data indicate the ability of the proposed method to reveal the non-continuity and fractures of surfaces with determining the dynamics size of the support domain based on the data structure. To compare the results of the proposed method with conventional interpolation methods including LPI, IDW, Kriging, and RBF, the root mean square error (RMSE), mean and median of errors are used. In this way, in addition to the overall accuracy of each method, the distribution of errors is also determined. The RMSE, mean and median errors of the proposed method, using the 10-fold cross-validation method for chromium (Cr), are 28.020, 0.2.201 and 2.874, respectively, and for iron (Fe) are 1.074, 0.017 and 0.094, respectively. Comparison of these results with conventional interpolation methods indicates the efficiency of the proposed method for both groups of high concentration and significant changes in the values and low concentration and almost uniform level of values. The results indicate the ability of the proposed method in detecting the jumps and non-continuity in the support domain and removal of some field points within the dynamic process, lead to a significant increase in the efficiency of the method compared to conventional methods.