We present results from numerical studies of neoclassical transport for multi‐species plasmas. The code, NEO, provides a first‐principles based calculation of the neoclassical transport coefficients directly from solution of the distribution function by solving a hierarchy of equations derived by expanding the fundamental drift‐kinetic equation in powers of ρ*i, the ratio of the ion gyroradius to system size. It extends previous studies by including the self‐consistent coupling of electrons and multiple ion species and strong toroidal rotation effects. Systematic calculations of the second‐order particle and energy fluxes and first‐order plasma flows and bootstrap current and comparisons with existing theories are given for multi‐species plasmas. The ambipolar relation ΣazaΓa = 0, which can only be maintained with complete cross‐species collisional coupling, is confirmed. The effects of plasma shaping are also explored.