The flux divergence technique of Athay and Skumanich (1967) is generalized for application to media whose properties vary in more than one spatial dimension. In this method, the flux divergence is viewed as an integro-differential functional of the source function. The source function is then expanded in terms of basis functions along characteristic paths, and, with the help of various interpolations, the flux divergence is converted to an approximate linear algebraic operator on a discrete spatial grid. A large but finite set of linear, inhomogeneous, simultaneous algebraic equations with known matrix coefficients is thus generated and is solved by direct matrix inversion for the source function at each point of the spatial grid.
Some aspects of the accuracy, stability, and computational convenience of the technique are discussed. Sample solutions for depth dependent, axially symmetric variations of temperature are shown.
