The finite difference (FD), finite element (FE), and finite volume (FV) methods are critically assessed by comparing the solutions produced by the three methods for a simple one-dimensional steady-state heat conduction problem with heat generation. Three issues are assessed: (1) accuracy of temperature, (2) accuracy of heat flux, and (3) satisfaction of global energy conservation. It is found that if the order of accuracy of the numerical discretization schemes is the same (central difference for FD and FV, linear basis functions for FE), the accuracy of the temperature produced by the three methods is similar, except close to the boundaries where the FV method outshines the other two methods. Consequently, the FV method is found to predict more accurate heat fluxes at the boundaries compared to the other two methods and is found to be the only method that guarantees both local and global conservation of energy irrespective of mesh size. The FD and FE methods both violate energy conservation, and the degree to which energy conservation is violated is found to be mesh size dependent. Furthermore, it is shown that in the case of prescribed heat flux (Neumann) and Newton cooling (Robin) boundary conditions, the accuracy of the FD method depends in large part on how the boundary condition is implemented. If the boundary condition and the governing equation are both satisfied at the boundary, the predicted temperatures are more accurate than in the case where only the boundary condition is satisfied.