Nonzero fluxes going through boundaries interfaces are normally observed in heat transfer, which in general can be described as inhomogeneous Neumann boundary conditions (BCs). Both smoothed particle hydrodynamics (SPH) and peridynamics have been employed for modeling heat transfer or thermal diffusion processes. The former is a numerical method used to approximate the solutions of classical heat diffusion PDEs. The latter provides a nonlocal model for heat diffusion. They both employ a nonlocal formulation, which requires a full support of the nonlocal kernel to ensure accuracy. In this work, we propose a new, higher-order method to enforce inhomogeneous Neumann BCs in SPH and peridynamic model for heat transfer problems. In that, fictitious layers of (ghost) particles are needed to guarantee full support of the nonlocal kernel. The temperature is extrapolated to the ghost particles based on the Taylor expansion and the BC to be imposed. By such, no additional term is introduced into the heat equation; meanwhile, the numerical solutions converge to the classical solutions with notably improved accuracy. To validate, assess, and demonstrate the proposed method, we simulate different transient or steady heat transfer problems subject to linear or nonlinear BCs, including heat conduction, natural convection, and presence of insulated cracks. The numerical results are compared with the exact solutions of classical PDEs, solutions of other numerical methods, or experimental data. (C) 2019 Elsevier Ltd. All rights reserved.