  We propose an efficient numerical algorithm to solve Bogoliubov de Gennes
equations self-consistently for inhomogeneous superconducting systems with a
reformulated polynomial expansion scheme. This proposed method is applied to
typical issues such as a vortex under randomly distributed impurities and a
normal conducting junction sandwiched between superconductors. With various
technical remarks, we show that its efficiency becomes remarkable in
large-scale parallel performance.
