Numerical simulations of binary neutron star mergers invariably show that, when a long-lived remnant forms, its rotation profile is never a simple decaying function of the radius but rather exhibits a maximum rotation rate shifted away from the center. This is in contrast to the usual differential rotation profile employed for the numerical modeling of axisymmetric equilibria of relativistic stars. Two families of rotation rate functions that mimic post-merger profiles were proposed by Uryū et al. (2017). In this work we implement Uryū’s profiles into the XNS code by Bucciantini and Del Zanna (2011) and we present novel equilibrium sequences of differentially rotating neutron stars. These are constructed by using three different equations of state, in order to study the dependence of mass, radius, angular momentum, and other important physical quantities, especially the quadrupole deformation and metric quadrupole moment, from the rotation properties.