Self-consistent approaches to superfluid many-fermion systems in three-dimensions (and their subsequent use in time-dependent studies) require a large number of diagonalizations of very large dimension Hermitian matrices, which results in enormous computational costs. We present an approach based on the shifted conjugate-orthogonal conjugate-gradient (COCG) Krylov method for the evaluation of the Green's function, from which we subsequently extract various densities (particle number, spin, current, kinetic energy, anomalous, etc.) of a nuclear system. The approach eschews the determination of the quasiparticle wavefunctions and their corresponding quasiparticle energies, which never explicitly appear in the construction of a single-particle Hamiltonian or needed for the calculation of various static nuclear properties, which depend only on densities. As benchmarks we present calculations for nuclei with axial symmetry, including the ground state of spherical (magic or semi-magic) and axially deformed nuclei, the saddle-point in the 240 Pu constrained fission path, and a vortex in the neutron star crust and demonstrate the superior efficiency of the shifted COCG Krylov method over traditional approaches.