Numerical solution of heat conduction in a heterogeneous material with small spatial and time scales can lead to excessive compute times due to the dense computational grids required. This problem is avoided by averaging the energy equation over the small-scales, which removes the appearance of the short spatial and time scales while retaining their effect on the average temperature. Averaging does, however, increase the complexity of the resulting thermal energy equation by introducing mixed spatial derivatives and six different averaged conductivity terms for three-dimensional analysis. There is a need for a numerical method that efficiently and accurately handles these complexities as well as the other details of the averaged thermal energy equation. That is the topic of this paper as it describes a numerical solution for the averaged thermal energy equation based on Fourier conduction reported recently in the literature. The solution, based on finite difference techniques that are second-order time-accurate and noniterative, is appropriate for three-dimensional time-dependent and steady-state analysis. Speed of solution is obtained by spatially factoring the scheme into an alternating direction sequence at each time level. Numerical stability is enhanced by implicit algorithms that make use of the properties of tightly banded matrices. While accurately accounting for the nonlinearity introduced into the energy equation by temperature-dependent properties, the numerical solution algorithm requires only the consideration of linear systems of algebraic equations in advancing the solution from one time level to the next. Computed examples are included and compared with those for a homogeneous material.