## 1 Introduction  
Since the 1940s, the knowledge of using composite laminate materials in industry has improved significantly [1]. A variety of methods have been developed to deal with composite laminate materials, which make the subject of composite laminate materials a matured discipline in applied science at present. The industrial use of composite laminate materials is widespread because they have the advantages of their anisotropic nature, which allows the material to be applied to a variety of engineering applications. However, composite materials may fail if subjected to severe environments such as high temperatures, even though there is no external load applied to the composites [2]. Specifically, failure in layered composites is caused by thermal induced stresses which are generated at the interface between different materials because the temperature distribution in the composites is non-uniform or discontinuous [3]. In order to understand how thermal stress is generated and distributed in the layered composites, it is important to investigate the temperature distribution in the composite laminate materials. While the mechanical behavior of composite laminate materials hasbeen explored for decades, studies on the heat conduction problems of layered composites using the meshless methods are less explored, which initiated this study. In the past, many heat conduction problems in layered composites have been routinely solved by numerical methods such as the finite difference method [4; 5; 6], finite element method [7; 8], finite volume method [9; 10; 11] and the boundary element method [12; 13]. In contrast to mesh-based numerical methods, the meshless methods, which do not need the mesh generation and boundary integral have been proposed such as the analytical method [14], method of fundamental solutions (MFS) [15; 16; 17], boundary knot method (BKM) [18], collocation Trefftz method (CTM) [19; 20; 21], radial basis function collocation method (RBFCM) [22; 23; 24; 25], element-free Galerkin method (EFG) [26], reproducing kernel particle method (RKPM) [27; 28], modified polynomial expansion method [29], meshless local boundary integral equation method (LBIE) [30; 31], and so on. Among these, boundary-type meshless methods have attracted considerable attention because of their simplicity. The CTM [32] is categorized as a boundary-type meshless method for the numerical solution of problems where approximate solutions are expressed as a truncated series of T-complete basis functions automatically satisfying governing equations. The use of the CTM is less widespread because the system of linear equations obtained from the Trefftz method is an ill-posed system [33]. On the other hand, the MFS is also a boundary-type meshless method for solving problems where the solutions are approximated by the fundamental solution which is expressed in terms of source points. The MFS requires placing source points outside the domain of the problem to avoid the effects of the singular characteristics of the fundamental solution; however, it often encounters difficulties such as finding an appropriate location for the source points [34; 35; 36]. In this study, we propose a novel hybrid boundary-type meshless method for solving two-dimensional stationary heat conduction problems in layered composite materials. The proposed method combines the collocation scheme from the MFS with the CTM to improve the applicability of both methods. Particular non-singular basis functions from the cylindrical harmonics are adopted in which the numerical solutions are approximated by superpositioning of the non-singular basis functions expressed in terms of many source points. For the modeling of multi-layer composite materials, we adopted the domain decomposition method (DDM) [37], which splits the domain into smaller subdomains that are intersected only at the interface between layers. For each subdomain, there exists an independent layer with its own thermal conductivity.