We propose a scheme based on a multidimensional Chebyshev interpolation to approximate smooth functions that depend on more than one variable. The present method generalizes the one dimensional Chebyshev approximation. The multidimensional approach can be used for generating databases like equation of state in the warm and hot dense matter. It is well suited to the present advance of massively parallel supercomputers.