## 3 Formulation of the Multiple Source Meshless Method  
We may express the heat equation in the dimensionless form as follows:  
\[\frac{\partial^{2}T}{\partial\rho^{2}}+\frac{1}{\rho}\frac{\partial^{2}T}{ \partial\rho^{2}}+\frac{1}{\rho^{2}}\frac{\partial^{2}T}{\partial\theta^{2}}=0 \tag{6}\]  
where \(\rho=\beta/R_{0}\). The radial coordinate is denoted by \(\rho\) and the angular coordinate is denoted by \(\theta\). Figure 3 displays the point distributions for the CTM, MFS, and MSMM in the simply-, doubly-, and multiply-connected regions. By using the additional theorem, approximation solutions can be obtained through the linear superposition principle of the particular solutions by arranging multiple sources within the domain, as shown in Figure 4. Then, the unknown coefficients \(a_{\rm jk}\) and \(b_{\rm jk}\) are determined using the boundary values as follows:  
\[T(\rho,\theta)\approx\sum\
olimits_{j=1}^{N_{p}}\sum\
olimits_{k=1}^{N_{q}}[a _{\rm jk}(\rho/R_{0})^{k}cosk\theta+b_{\rm jk}\Big{(}\rho/R_{0})^{k}sink\theta \Big{)}, \tag{7}\]  
\[\frac{\partial T(\rho,\theta)}{\partial n}\approx\sum\
olimits_{j=1}^{N_{p}} \sum\
olimits_{k=1}^{N_{q}}[a_{\rm jk}\frac{\partial\big{(}(\rho/R_{0})^{k} cosk\theta\big{)}}{\partial n}+b_{\rm jk}\frac{\partial\big{(}(\rho/R_{0})^{k}sink \theta\big{)}}{\partial n}], \tag{8}\]  
where \(a_{\rm jk}\) and \(b_{\rm jk}\) denote unknown coefficients, \(N_{p}\) denotes the source point number, and \(N_{q}\) denotes the number of Trefftz terms. The characteristic length \(R_{0}\) is adopted from [31] as:  
\[R_{0}=\xi\times max(\rho), \tag{9}\]  
where \(\xi\) denotes a constant that is a user-defined parameter and \(max(\rho)\) represents the maximum distance between the boundary and source points of the problem. By substituting Equations (7) and (8) into Equations (2) and (3), we obtain the Dirichlet and Neumann boundary conditions as follows:  
\[T(\rho,\theta)\approx\sum\
olimits_{j=1}^{N_{p}}\sum\
olimits_{k=1}^{N_{q}} \Big{[}a_{jk}C_{k}+b_{jk}D_{k}\Big{]}, \tag{10}\]  
\[\frac{\partial T(\rho,\theta)}{\partial n}\approx\sum\
olimits_{j=1}^{N_{p}} \sum\
olimits_{k=1}^{N_{q}}\Big{[}a_{jk}E_{k}+b_{jk}F_{k}\Big{]}, \tag{11}\]  
where \(C_{k}:=(\frac{\rho}{R_{0}})^{k}cosk\theta,D_{k}:=(\frac{\rho}{R_{0}})^{k}sink \theta,E_{k}:=\frac{\partial((\frac{\rho}{R_{0}})^{k}\cosk\theta)}{\partial n}\), and \(F_{k}:=\frac{\partial((\frac{\rho}{R_{0}})^{k}sink\theta)}{\partial n}\). By using Equations (10) and (11), we obtain a linear equation system as follows:  
\[A\alpha=b. \tag{12}\]  
Then, Equation (12) can be rewritten as follows:  
\[A=\left[\begin{array}{cccc}D_{1,k=1}^{l}&D_{1,k=1}^{ll}&\ldots&D_{1,k=Nq}^{ l}&D_{1,k=Nq}^{ll}\\ D_{2,k=1}^{l}&D_{2,k=1}^{ll}&\ldots&D_{2,k=Nq}^{l}&D_{2,k=Nq}^{ll}\\ D_{3,k=1}^{l}&D_{3,k=1}^{ll}&\ldots&D_{3,k=Nq}^{l}&D_{3,k=Nq}^{ll}\\ \vdots&\vdots&\ldots&\vdots&\vdots\\ D_{i,k=1}^{l}&D_{j,k=1}^{ll}&\ldots&D_{i,k=Nq}^{l}&D_{i,k=Nq}^{ll}\\ N_{1,k=1}^{ll}&N_{1,k=1}^{ll}&\ldots&N_{1,k=Nq}^{l}&N_{1,k=Nq}^{ll}\\ N_{2,k=\[N_{jk}^{l}:= \begin{array}{l}[(k(1/R_{0})^{k}\rho_{j}^{k-1}cos(k\theta_{j})cos\theta_{ j})+(k(1/R_{0})^{k}\rho_{j}^{k-1}sin(k\theta_{j})sin\theta_{j})]n_{x}\\ +[(k(1/R_{0})^{k}\rho_{j}^{k-1}cos(k\theta_{j})sin\theta_{j})-(k(1/R_{0})^{k} \rho_{j}^{k-1}sin\left(k\theta_{j}\right)cos\theta_{j})]n_{y}\\ N_{jk}^{ll}:= \begin{array}{l}[(k(1/R_{0})^{k}\rho_{j}^{k-1}sin(k\theta_{j})cos\theta_{ j})-(k(1/R_{0})^{k}\rho_{j}^{k-1}cos(k\theta_{j})sin\theta_{j})]n_{x}\\ +[(k\left(1/R_{0})^{k}\rho_{j}^{k-1}sin\left(k\theta_{j}\right)sin\theta_{j} \right)+\left(k\left(1/R_{0}\right)^{k}\rho_{j}^{k-1}cos\left(k\theta_{j} \right)cos\theta_{j}\right)]n_{y}\\ \end{array}.\end{array} \]  
Here, \(A\) denotes the coefficient matrix, \(\alpha\) denotes the unknown coefficients and is given as \(\alpha=[a_{0},\dots,\ N_{q}]^{T}\), and \(b\) denotes the function value from the boundary conditions. The DDM [32] is adopted for problems in the doubly-connected domain.