Even in the Newtonian approximation, analytic models of rotating stars with specified equations of state are limited to the classical uniform-density ellipsoids, and analytic relativistic models are restricted to rotating disks [469]. But in the full theory, as in the Newtonian framework, accurate numerical models are readily obtained. All current codes involve an iteration in which one solves a set of nonlinear equations for the metric functions, together with the equation of hydrostationary equilibrium. In nearly all codes, one successively solves the equations for the metric functions for one function at a time, with the others set to their values from the last iteration; and the first integral (1.110) of the equation of hydrostationary equilibrium (with the updated metric functions) is then used to update the matter distribution. The methods fall into two classes, depending on whether each nonlinear equation is solved by a Newton-Raphson method or, in the case of elliptic equations for the metric functions, by inverting a flat-space elliptic operator with known Green's function. For uniform rotation and a specified barotropic equation of state, the family of equilibria is 2-dimensional (locally specified, for example, by central density and angular velocity Ω), and within each class of code, different choices are made of the two parameters that are held fixed to obtain a convergent iteration. Finally, the codes differ in whether spectral methods are used to represent radial and angular functions and to evaluate corresponding integrals and derivatives, and in the choice of independent components of the field equation that are to be solved.