  We have developed a theoretical model and a numerical code for stationary
rotating superfluid neutron stars in full general relativity. The underlying
two-fluid model is based on Carter's covariant multi-fluid hydrodynamic
formalism. The two fluids, representing the superfluid neutrons on one hand,
and the protons and electrons on the other, are restricted to uniform rotation
around a common axis, but are allowed to have different rotation rates. We have
performed extensive tests of the numerical code, including quantitative
comparisons to previous approximative results for these models. The results
presented here are the first ``exact'' calculations of such models in the sense
that no approximations (other than that inherent in a discretized numerical
treatment) are used. Using this code we reconfirm the existence of
prolate-oblate shaped configurations. We studied the dependency of the Kepler
rotation limit and of the mass-density relation on the relative rotation rate.
We further demonstrate how one can simulate a (albeit fluid) neutron-star
``crust'' by letting one fluid extend further outwards than the other, which
results in interesting cases where the Kepler limit is actually determined by
the outermost but slower fluid.
