Abstract The method of fundamental solutions (MFS) belongs to the family of meshless boundary collocation methods and now has been successfully tried for many kinds of engineering problems. The traditional MFS based on the “global” boundary discretization, however, leads to dense and non-symmetric coefficient matrices that, although smaller in sizes, require huge computational cost to compute the system of equations using direct solvers. Such an approach will be arduous, time consuming and computationally expensive for analyzing large-scale problems. In the present work, a localized version of the MFS, named as the localized MFS (LMFS), is proposed for large-scale modelling of three-dimensional (3D) anisotropic heat conduction problems. In the LMFS, the computational domain can be divided into small subdomains with a simple geometry such as circle and/or sphere. To each of the subdomains, the MFS formulation is applied and the unknown coefficients on the local simple geometric boundary are approximated by the moving least square (MLS) method. The satisfactions of governing equations at interior points and boundary conditions at boundary nodes lead to a sparse and banded system matrix. Numerical examples with up to 1,000,000 unknowns are solved successfully using the developed LMFS code. The advantages, disadvantages and potential applications of the proposed method, as compared with the traditional MFS and boundary element method (BEM), are discussed. Finally, a fast, reliable and self-contained MATLAB code is provided in the part of Supplementary Materials of the paper.