[en] It is possible to solve three-dimensional steady state thermoelastic problems from data of temperature, fluxes, tractions, and displacements on the boundary alone. The numerical discretisation of these integral equations is based on a technique similar to that of finite elements. The boundary alone needs to be discretised. Two-dimensional examples are given to prove the accuracy and the feasibility. In the case of three-dimensional problems, the surface is represented by eight nodes quadrilateral elements and six nodes triangular elements. The unknown (temperature, flux, displacement, traction) may be considered to vary linearly, quadratically, or cubically, with respect to he intrinsic coordinates of each element. The integration is performed numerically using Gaussian quadrature formulas, for which the number of integration points is chosen automatically by the program, so that the upper bound of error in integration is minimized. In order to obtain a banded form matrix and also to be able to study elongated structures, the body is divided into subregions, for each of which the integral equations are written. The thermoelastic stress and deformation field are obtained from two successive calculations on the same mesh. The first gives the thermal field. The second, taking account of the thermal field, calculates the thermoelastic displacement and stress field. This type of approach is especially suited to complicated three-dimensional thick structures for which finite element procedures are very expensive. The data are simple to generate, for data on the boundary alone have to be given. Much time can be saved by the user of the program in data generations and data checks