In the paper, we point out a drawback of the Fourier sine series method to represent a given odd function, where the boundary Gibbs phenomena would occur when the boundary values of the function are non-zero. We modify the Fourier sine series method by considering the consistent conditions on the boundaries, which can improve the accuracy near the boundaries. The modifications are extended to the Fourier cosine series and the Fourier series. Then, novel boundary consistent methods are developed to solve the 1D and 2D heat equations. Numerical examples confirm the accuracy of the boundary consistent methods, accounting for the non-zeros of the source terms and considering the consistency of heat equations on the boundaries, which can not only overcome the near boundary errors but also improve the accuracy of solution about four orders in the entire domain, upon comparing to the conventional Fourier sine series method and Duhamel’s principle.