This paper presents a systematic approach to determine temperature wall functions for high Rayleigh number flows using asymptotics. An asymptotic analysis of the flow and heat transfer in the near wall region forms the basis for the development of the wall functions. Appropriate normalization of the variables followed by asymptotic matching of the temperature gradients of the inner and outer layers in the overlap region leads to a logarithmic temperature profile as a wall function that has undetermined constants. A key classification that has been made in the present study is the introduction of (1) The direct problem and (2) The inverse problem. The former means that temperature profiles, either from experiments or Direct Numerical Simulations (DNS), are available and the wall function problem finally reduces to one of determining certain constants in a general wall function formula. More radical and of more interest, is the inverse problem. The idea behind this it is that when a temperature profile can be recast into a Nusselt-Rayleigh correlation, it should be perfectly possible for one to start from a Nusselt-Rayleigh correlation and end up with a wall function for temperature. This approach again will have undetermined constants that can be calibrated from either experimental or DNS data. The main advantage of using the inverse problem is the dispensation of the need to measure temperatures accurately within the boundary layer. For both the direct and inverse problems, a graded treatment to determine the constants is presented. The treatment at its highest level will result in a parameter estimation problem that can be posed as an optimization problem. The optimization problem is then solved by state of the art techniques like Levenberg-Marquardt algorithm and Genetic algorithms (GA) and the solutions are compared. While for the direct problem, the approach is illustrated for the infinite channel problem (a simple flow), for the inverse problem, the approach is elucidated for the Rayleigh-Benard problem (a complex flow). Finally, a blending procedure to arrive at a universal temperature profile that is valid in the viscous sublayer, buffer and the overlap layers is suggested. The key ideas of (1) using optimization techniques for determining the constants in the wall function and (2) obtaining wall functions from the Nusselt numbers by the inverse approach are expected to be useful for a wide class of problems involving natural convection. (c) 2007 Elsevier Ltd. All rights reserved.