Heat diffusion has been widely used in brain imaging for surface fairing, mesh regularization and noisy cortical data smoothing. In the previous spectral decomposition of graph Laplacian, Chebyshev polynomials were only used. In this paper, we present a new general spectral theory for the Laplace-Beltrami operator on a manifold that works for an arbitrary orthogonal polynomial with a recurrence relation. Besides the Chebyshev polynomials that was previous used in diffusion wavelets and convolutional neural networks, we provide three other polynomials to show the generality of the method. We also derive the closed-form solutions to the expansion coefficients of the spectral decomposition of the Laplace-Beltrami operator and use it to solve heat diffusion on a manifold for the first time. The proposed fast polynomial approximation scheme avoids solving for the eigenfunctions of the Laplace-Beltrami operator, which are computationally costly for large mesh size, and the numerical instability associated with the finite element method based diffusion solvers. The proposed method is applied in localizing the male and female differences in cortical sulcal and gyral graph patterns obtained from MRI.