An Eulerian numerical discretization scheme for the solution of the first-order drift-kinetic equation is presented in detail. The approach is valid for multi-species plasmas, including complete impurity and electron physics, and thus ambipolarity is properly maintained. The code, NEO, provides a complete description of the second-order neoclassical transport fluxes and first-order flows, including the effects of strong (finite Mach number) toroidal rotation. Corrections to the weak rotation limit are demonstrated for multi-species plasmas over a wide range of collisionality.