Interpolation of the function of two variables with large gradients in the boundary layers

The problem of interpolation for the function of two variables with large gradients in the boundary layers is investigated. It is suggested that the function can be represented as a sum of a regular component with bounded derivatives up to some order and of two boundary layer components. The boundary layer components are known, but their coefficients are uncertain. Such representation is typical for the solution of a singular perturbed elliptic problem. A two-dimensional interpolation formula, which is exact on the boundary layer components, is deduced. The formula has the arbitrary number of nodes in each direction. The accuracy estimate, which is uniform in gradients of the interpolated function in the boundary layers, is proved. The results of numerical experiments are provided.