Partial differential equations arising in engineering and other sciences describe nature adequately in terms of symmetry properties. This article develops a numerical method based on the Laplace transform and the numerical inverse Laplace transform for numerical modeling of diffusion problems. This method transforms the time-dependent problem to a corresponding time-independent inhomogeneous problem by employing the Laplace transform. Then a local radial basis functions method is employed to solve the transformed problem in the Laplace domain. The main feature of the local radial basis functions method is the collocation on overlapping sub-domains of influence instead of on the whole domain, which reduces the size of the collocation matrix; hence, the problem of ill-conditioning in global radial basis functions is resolved. The Laplace transform is used in comparison with a finite difference technique to deal with the time derivative and avoid the effect of the time step on numerical stability and accuracy. However, using the Laplace transform sometimes leads to a solution in the Laplace domain that cannot be converted back into the real domain by analytic methods. Therefore, in such a case, the Laplace transform is inverted numerically. In this investigation, two inversion techniques are utilized; (i) the contour integration method, and (ii) the Stehfest method. Three test problems are used to evaluate the proposed numerical method. The numerical results demonstrate that the proposed method is computationally efficient and highly accurate.