We demonstrate a method to systematically obtain eigenvalues and eigenstates of a many-body Hamiltonian describing collective neutrino oscillations. The method is derived from the Richardson-Gaudin framework, which involves casting the eigenproblem as a set of coupled nonlinear “Bethe Ansatz equations”, the solutions of which can then be used to parametrize the eigenvalues and eigenvectors. The specific approach outlined in this paper consists of defining auxiliary variables that are related to the Bethe-Ansatz parameters, thereby transforming the Bethe-Ansatz equations into a different set of equations that are numerically better behaved and more tractable. We show that it is possible to express not only the eigenvalues, but also the eigenstates, directly in terms of these auxiliary variables without involving the Bethe Ansatz parameters themselves. In this paper, we limit ourselves to a two-flavor, single-angle neutrino system.