We present a new discrete chemo-dynamical axisymmetric modeling technique, which we apply to the dwarf spheroidal galaxy Sculptor. The major improvement over previous Jeans models is that realistic chemical distributions are included directly in the dynamical modelling of the discrete data. This avoids loss of information due to spatial binning and eliminates the need for hard cuts to remove contaminants and to separate stars based on their chemical properties. Using a combined likelihood in position, metallicity and kinematics, we find that our models naturally separate Sculptor stars into a metal-rich and a metal-poor population. Allowing for non-spherical symmetry, our approach provides a central slope of the dark matter density of $gamma = 0.5 pm 0.3$. The metal-rich population is nearly isotropic (with $beta_r^{red} = 0.0pm0.1$) while the metal-poor population is tangentially anisotropic (with $beta_r^{blue} = -0.2pm0.1$) around the half light radius of $0.26$ kpc. A weak internal rotation of the metal-rich population is revealed with $v_{max}/sigma_0 = 0.15 pm 0.15$. We run tests using mock data to show that a discrete dataset with $sim 6000$ stars is required to distinguish between a core ($gamma = 0$) and cusp ($gamma = 1$), and to constrain the possible internal rotation to better than $1,sigma$ confidence with our model. We conclude that our discrete chemo-dynamical modelling technique provides a flexible and powerful tool to robustly constrain the internal dynamics of multiple populations, and the total mass distribution in a stellar system.