In recent years, radial basis function (RBF) based meshless methods have been used extensively to solve partial differential equations (PDEs) arising in thermal sciences. Meshless methods (RBF) make use of multi-variable interpolation on scattered set of data points to approximate a variable at a defined location and do not need grid connectivity information in the form of faces or control volumes as required by grid based techniques such as finite difference (FDM), finite volume (FVM), boundary element (BEM) methods, etc. Over the years, it has been demonstrated that polyharmonic splines (PHS-RBF) with appended polynomial can be used as an alternative to the conventional radial basis functions (RBFs) as PHS-RBF is independent of shape parameter and the accuracy or convergence can be controlled to a very large extent by the appended polynomial. In the present work, we demonstrate the capability of PHSRBF to solve heat conduction problems in multiple twodimensional geometries with emphasis on the discretization error. Along with the discretization error, a local condition number analysis is also presented in detail. An improved accuracy and convergence rate is achieved with the increase in the number of nodes and degree of appended polynomial. PHS-RBF has proven to be a robust numerical technique to solve heat conduction problems in two-dimensional (2D) geometries.