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.