Summary: The article discusses the implementation of the algorithm for solving axisymmetric contact problems of the thermoelasticity theory using the mortar method. This algorithm is used for numerical simulation of the contact interaction of several bodies under thermomechanical loading. The ill-conditioned system of linear algebraic equations obtained as a result of finite element discretization is numerically solved using the modified symmetric successive overrelaxation method (MSSOR), generalized to the case of contact of several bodies. The results of the algorithm application are demonstrated on a problem simulating some processes in a fuel element with a different number of bodies. The effect of the contacting bodies number and mesh steps on the number of iterations necessary to achieve a given accuracy while solving the system of equations is investigated.