This paper focuses on implementation of the sampling surfaces (SaS) method for the three-dimensional (3D) thermal stress analysis of steady-state thermoelasticity problems for laminated functionally graded (FG) shells. The SaS formulation is based on choosing inside the nth layer I n not equally spaced SaS parallel to the middle surface of the shell in order to introduce the temperatures and displacements of these surfaces as basic shell variables. Such choice of unknowns permits the presentation of the proposed thermoelastic FG shell formulation in a very compact form. The SaS are located inside each layer at Chebyshev polynomial nodes that improves the convergence of the SaS method significantly. As a result, the SaS formulation can be applied efficiently to analytical solutions for laminated FG shells, which asymptotically approach the 3D exact solutions of thermoelasticity as the number of SaS I n tends to infinity.