We report an investigation of temperature profiles in turbulent Rayleigh–Bénard convection of water based on direct numerical simulations (DNS) for a cylindrical cell with unit aspect ratio for the same Prandtl number Pr and similar Rayleigh numbers Ra as used in recent high-precision measurements by Funfschilling et al. (J. Fluid Mech., vol. 536, 2005, p. 145). The Nusselt numbers Nu computed for Pr = 4.38 and Ra = 108, 3 × 108, 5 × 108, 8 × 108 and 109 are found to be in excellent agreement with the experimental data corrected for finite thermal conductivity of the walls. Based on this successful validation of the numerical approach, the DNS data are used to extract vertical profiles of the mean temperature. We find that near the heating and cooling plates the non-dimensional temperature profiles Θ(y) (where y is the non-dimensional vertical coordinate), obey neither a logarithmic nor a power law. Moreover, we demonstrate that the Prandtl–Blasius boundary layer theory cannot predict the shape of the temperature profile with an error less than 7.9% within the thermal boundary layers (TBLs). We further show that the profiles can be approximated by a universal stretched exponential of the form Θ(y) ≈ 1 − exp(−y − 0.5y2) with an absolute error less than 1.1% within the TBLs and 5.5% in the whole Rayleigh cell. Finally, we provide more accurate analytical approximations of the profiles involving higher order polynomials in the approximation.