There has been increasing interest in studying the Richardson model from which one can derive the exact solution for certain pairing Hamiltonians. However, it is still a numerical challenge to solve the nonlinear equations involved. In this paper we tackle this problem by employing a simple hybrid polynomial approach. The method is found to be robust and is valid for both deformed and nearly spherical nuclei. It also provides important and convenient initial guesses for spherical systems with large degeneracy. As an example, we apply the method to study the shape coexistence in neutron-rich Ni isotopes.