  We study models of uniformly and differentially rotating neutron stars in the
framework of post-Newtonian approximation in general relativity as established
by Chandrasechar. In particular, we adopt the polytropic equation of state in
order to derive the appropriate hydrodynamic equations and a rotation law based
on the generalized Clement's model. To compute equilibrium configurations at
the mass-shedding limit, i.e. at critical angular velocity (equivalently,
Keplerian angular velocity), we develop an iterative numerical method,
belonging to the category of the well-known ``self-consistent field methods'',
with two perturbation parameters: the ``rotation parameter'' $\bar{\upsilon}$
and the ``gravitation or relativity parameter'' $\bar{\sigma}$. These two
parameters represent the effects of rotation and gravity on the configuration.
We investigate the validity and the limits of our method by comparing our
results with respective results of other computational methods and public
domain codes. As it turns out, our method can derive satisfactory results for
general-relativistic polytropic configurations at critical rotation.
