A node-based smoothed point interpolation method (NS-PIM) is formulated to analyze steady-state thermoelastic problems. In this approach, shape functions are constructed using the point interpolation method (PIM). which permits the straightforward enforcement of essential boundary conditions. The smoothed Galerkin weak form is then used to construct discretized system equations using smoothing domains constructed based on nodes. The bound property, accuracy and convergence of the present formulation are studied using 1D and 2D thermoelasticity problems. It is found that the computed temperature and its resulted gradient are in very good agreement with the analytical results or those obtained using the finite element method (FEM). Compared with the 3-node triangular FEM. the NS-PIM can achieve better accuracy and higher convergence in energy norm using the same linear triangular mesh. Together with the FEM, we now for the first time have a simple way to obtain both upper and lower bounds of the exact solution to thermoelasticity problems. Crown Copyright (C) 2008 Published by Elsevier Ltd. All rights reserved.