The three-dimensional problem of the theory of elasticity related to isotropic layer compression by normal load, distributed on a limited area, is solved by the method of initial functions (MIF). The layer divided into separate sub layers with different elastic characteristics serves as a model of the multilayer foundation. A parallelepiped cut out from the infinite layer with dimensions much larger than those of the load area may be considered as an elastic half-space. A numerical-analytic solution was obtained by a specially designed program based on the symbolic computation system called Maple, in which the desired functions of the displacements are represented by a Fourier series. The problem of computational instability calculations inherent in MIF at high numbers of harmonics was solved by using representation of real numbers with sufficient length mantissa. The results were compared both with solutions of the classical theory of elasticity for elastic halfspace, stipulated in the guidelines for designing foundations, and with the finite element method solutions.