This paper is to present the treatment of internal thermoelastic stress analysis in 3D anisotropic bodies by the boundary element method (BEM). Fundamentally, thermal effects will give rise to an additional volume integral in the boundary integral equation (BIE). By applying the fundamental solutions represented by Fourier series, the volume integral has been analytically transformed to the boundary. For the present work, spatial differentiations of the integral equation are performed to give displacement gradients at internal points of interest. This differentiated integral equation is further implemented to perform thermoelastic stress analysis inside 3D anisotropic bodies. This analysis is particularly important in engineering applications when thermoelastic stresses concentrations are present inside the bodies. The present work is the first BEM implementation for this study by the transformed BIE. In the end, two benchmark examples are tested to demonstrate the applicability of the present BEM treatment.