  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.
