  We propose a highly efficient numerical method to describe inhomogeneous
superconductivity by using the kernel polynomial method in order to calculate
the Green's functions of a superconductor. Broken translational invariance of
any type (impurities, surfaces or magnetic fields) can be easily incorporated.
We show that limitations due to system size can be easily circumvented and
therefore this method opens the way for the study of scenarios and/or
geometries that were unaccessible before. The proposed method is highly
efficient and amenable to large scale parallel computation. Although we only
use it in the context of superconductivity, it is applicable to other
inhomogeneous mean-field theories.
