A new spectral Jacobi–Gauss–Lobatto collocation (J–GL–C) method is developed and analyzed to solve numerically parabolic partial differential equations (PPDEs) subject to initial and nonlocal boundary conditions. The method depends basically on the fact that an expansion in a series of Jacobi polynomials is assumed, for the function and its space derivatives occurring in the partial differential equation (PDE), the expansion coefficients are then determined by reducing the PDE with its boundary conditions into a system of ordinary differential equations (SODEs) for these coefficients. This system may be solved numerically in a step-by-step manner by using implicit the Runge–Kutta (IRK) method of order four. The proposed method, in contrast to common finite-difference and finite-element methods, has the exponential rate of convergence for the spatial discretizations. Numerical results indicating the high accuracy and effectiveness of this algorithm are presented.