  A highly accurate, multi-domain spectral code is used in order to construct
sequences of general relativistic, differentially rotating neutron stars in
axisymmetry and stationarity. For bodies with a spheroidal topology and a
homogeneous or an N=1 polytropic equation of state, we investigate the solution
space corresponding to broad ranges of degree of differential rotation and
stellar densities. In particular, starting from static and spherical
configurations, we analyse the changes of the corresponding surface shapes as
the rate of rotation is increased. For a sufficiently weak degree of
differential rotation, the sequences terminate at a mass-shedding limit, while
for moderate and strong rates of differential rotation, they exhibit a
continuous parametric transition to a regime of toroidal fluid bodies. In this
article, we concentrate on the appearance of this transition, analyse in detail
its occurrence and show its relevance for the calculation of astrophysical
sequences. Moreover, we find that the solution space contains various types of
spheroidal configurations, which were not considered in previous work, mainly
due to numerical limitations.
