The paper deals with inverse thermoelasticity problems. A general, robust, and numerically efficient technique, for retrieving the multi-dimensional, highly varying distributions of boundary conditions is presented. In the class of considered inverse problems, both input data and sought-for quantities are usually specified at the domain boundary, only. As a sequence of forward sub-solutions underlies the inverse analysis, the numerical method of choice for solving the field problem is the boundary element method (BEM). The derived inverse technique is capable of retrieving boundary condition distributions in steady state and transient problems. The accuracy and stability of the algorithm are verified by considering problems involving constant, functionally graded, and temperature dependent material properties. Strain components and temperatures, subject to uncertainties, are used as input data. Presented numerical examples show that the method is capable of reconstructing mechanical and thermal loads with reasonable accuracy.