Axisymmetric steady solutions of Taylor-Couette flow, known as Taylor vortices (Taylor 1923; Davey 1962), are studied numerically and theoretically focusing on the high Taylor number regime. As the axial period of the solution shortens from about the gap, the Nusselt number (the torque on the cylinder walls normalised by its laminar value) goes through two peaks before returning to laminar flow. In this process, the asymptotic nature of the solution changes in four stages, as revealed by the asymptotic analysis. When the aspect ratio of the roll cell is about unity, the solution quantitatively well captures the characteristics of the classical turbulence regime in the experiments and DNSs summarised in Grossmann et al. (2016). Theoretically, the Nusselt number of the solution is proportional to the quarter power of the Taylor number. It is surprising that a single unstable steady-state solution can approximate turbulence to some extent. On the other hand, the maximised Nusselt number obtained by shortening the axial period can reach the experimental value around the onset of the ultimate turbulence regime, although at higher Taylor numbers the theoretical predictions eventually underestimate the experimental values. This suggests that ultimate turbulence has a mechanism that cannot be represented by steady axisymmetric flow. An important consequence of the asymptotic analyses is that the mean angular momentum should become uniform in the core region unless the axial period is too short. The theoretical scaling laws deduced for the steady solutions can be carried over to Rayleigh-B\'enard convection.