Advanced technological and engineering solutions involve usage of complex tribological systems. There is a demand for precise and computationally efficient methods to describe such systems. The aim of this study is to develop a computationally efficient method to solve the problem of deformation of a two-layer system using the Green's function molecular dynamics (GFMD) technique. We consider a viscoelastic layer attached to an elastic halfspace and derive a constitutive equation in Fourier space from the corresponding elastic solution. This third-order equation is numerically integrated by the backward Euler method, and a quasi-static solution is found through the fast inertial relaxation engine (FIRE) optimization algorithm. The method is illustrated with a simple model of indentation by a rigid cylinder. Using this method, contact area and pressure were calculated as a function of time for various shear modulus values.