  We present an efficient algorithm for obtaining the gauge-invariant gradient
expansion of the local density of states and the free energy of a clean
superconductor. Our method is based on a new mapping of the semiclassical
linearized Gorkov equations onto a pseudo-Schroedinger equation for a
three-component wave-function psi(x), where one component is directly related
to the local density of states. Because psi(x) satisfies a linear equation of
motion, successive terms in the gradient expansion can be obtained by simple
linear iteration. Our method works equally well for real and complex order
parameter, and in the presence of arbitrary external fields. We confirm a
recent calculation of the fourth order correction to the free energy by
Kosztin, Kos, Stone and Leggett [Phys. Rev. B 58, 9365 (1998)], who obtained a
discrepancy with an earlier result by Tewordt [Z. Phys. 180, 385 (1964)]. We
also give the fourth order correction to the local density of states, which has
not been published before.
