There are many different types of numerical analysis for partial differential equations. Numerical methods that have been studied over a long period of time and have been used in engineering practice include the finite difference method (FDM), finite element method (FEM), and boundary element method (BEM). This paper focuses on the FDM. Conventionally, numerical calculations using the FDM have mostly been performed with second-order accuracy. It usually guarantees a calculation accuracy of 3 to 4 significant figures, which is considered to be sufficient accuracy from an engineering perspective. We generally recognize that numerical calculations are approximate calculations. However, when the FDM is generalized in a way that incorporates the conventional theory, it becomes possible to perform unlimitedly high-accuracy numerical calculations, and it becomes possible to easily perform calculations that converge to a theoretical numerical solution with a certain number of significant figures. This type of numerical calculation system is defined as the interpolation FDM (IFDM). The essential requirement for the establishment of IFDM is the algebraic polynomial approximation of the (real number) analytic function. This paper describes this concept and outlines the theoretical features and computational examples of IFDM, which has three characteristics: (i) the ability to handle arbitrary domains and arbitrary boundary conditions, (ii) unlimited high-accuracy calculations, and (iii) high-speed calculations.