A numerical procedure is disclosed to improve the prediction of heat fronts when simulating hot fluid injection in viscous hydrocarbon reservoirs. The mathematical model is composed of the conventional governing equations that describe multiphase fluid flow and energy balance. The reservoir geometry can be partitioned into a regular Cartesian grid or an irregular corner-point geometry grid. The numerical procedure uses the finite different (FD) method to solve the flow equations and the discontinuous Galerkin (DG) method to solve the energy balance equation. The proposed FD-DG method is an alternative to the traditional solution procedure that uses the FD method to solve both the flow and the energy equations. The traditional method has the deficiency that it may require excessive number of grid cells to achieve acceptable resolution of the heat fronts. The proposed FD-DG method significantly reduces numerical dispersion near discontinuities in the solution of the energy equation and therefore provides a better capture of the heat fronts. To obtain a desired accuracy in the energy equation solution, the FD-DG method can be orders of magnitude faster than the traditional method. The superiority of the FD-DG method is that it converges on coarser grids while the traditional method requires much finer grids.