  We show how multi-level BCS Hamiltonians of finite systems in the strong
pairing interaction regime can be accurately approximated using
multi-dimensional shifted harmonic oscillator Hamiltonians. In the Shifted
Harmonic Approximation (SHA), discrete quantum state variables are approximated
as continuous ones and algebraic Hamiltonians are replaced by differential
operators. Using the SHA, the results of the BCS theory, such as the gap
equations, can be easily derived without the BCS approximation. In addition,
the SHA preserves the symmetries associated with the BCS Hamiltonians. Lastly,
for all interaction strengths, the SHA can be used to identify the most
important basis states -- allowing accurate computation of low-lying
eigenstates by diagonalizing BCS Hamiltonians in small subspaces of what may
otherwise be vastly larger Hilbert spaces.
