Multivariate interpolation is a fundamental and long-studied problem, which has numerous applications in mathematics, engineering, computer science, and the natural sciences. A basic tool for solving the problem of high-dimensional interpolation is through the usage of radial basis functions (RBFs). In fact, the combination of interpolation and RBFs can lead to very good properties in high dimensions. Unfortunately, the linear system of equations derived from the approximations of RBFs with a high order of convergence is ill-conditioned and unstable, and usually includes a full interpolation matrix. To solve such a system of equations, we face a very high computational cost if the dimension of the problems or the number of data points is increased. This can also lead to intense instability in the considered problem, and the condition number of the system of equations, as a measure of the ill-conditioning criterion, will be very large. To overcome these problems, this paper presents a layer-by-layer interpolation approach for solving scattered data approximation of an unknown multivariate function, where the information of this approximation problem has been given in certain layers. In this approach, by creating a layered structure, the computational cost is reduced and decreasing the condition number is also possible. This structure provides a possibility for encountering a much smaller linear system of equations. In other words, it increases the accuracy and the stability of the numerical structure of the considered interpolation. The new method can ensure the existence and the uniqueness of the solution for the multilayer interpolation problem. We also find that the layer-by-layer approach provides more numerically stable calculations than the traditional interpolation method. The new approach is applied for certain numerical examples in high dimensions and the obtained results confirm the high accuracy and the low computational cost of the proposed method. Finall