Simulations of homogeneously sheared turbulence (HST) are conducted until a universal self-similar state is established at the long non-dimensional time $\unicode[STIX]{x1D6E4}t=20$, where $\unicode[STIX]{x1D6E4}$ is the shear rate. The simulations are enabled by a new robust and discretely conservative algorithm. The method solves the governing equations in physical space using the so-called shear-periodic boundary conditions. Convection by the mean homogeneous shear flow is treated implicitly in a split step approach. An iterative Crank–Nicolson time integrator is chosen for robustness and stability. The numerical strategy captures without distortion the Kelvin modes, rotating waves that are fundamental to homogeneously sheared flows and are at the core of rapid distortion theory. Three direct numerical simulations of HST with the initial Taylor scale Reynolds number $Re_{\unicode[STIX]{x1D706}0}=29$ and shear numbers of $S_{0}^{\ast }=\unicode[STIX]{x1D6E4}q^{2}/\unicode[STIX]{x1D716}=3$, 15 and 27 are performed on a $2048\times 1024\times 1024$ grid. Here, $\unicode[STIX]{x1D716}$ is the dissipation rate and $1/2q^{2}$ is the turbulent kinetic energy. The long integration time considered allows the establishment of a self-similar state observed in experiments but often absent from simulations conducted over shorter times. The asymptotic state appears to be universal with a long time production to dissipation rate ${\mathcal{P}}/\unicode[STIX]{x1D716}\sim 1.5$ and shear number $S^{\ast }\sim 10$ in agreement with experiments. While the small scales exhibit strong anisotropy increasing with initial shear number, the skewness of the transverse velocity derivative decreases with increasing Reynolds number.