In case of natural convection modeling, when Boussinesq assumption is used, we encounter coupled nonlinear differential equations. In this work, the authors have modeled natural heat convection by implementing one of the newest analytical methods of solving nonlinear differential equations called homotopy analysis method (HAM), which gives us a vast freedom to choose the answer type. We have used an iterating analytical method in order that cope with nonlinearity. Also, we apply some provisions because of particular difficulties that are caused by coupling problem. A new adapting boundary condition is proposed in this work that is based on an initial guess and then it is developed to the solution expression. We must notice that HAM is applied to our case study according to the physics of the target problem. (C) 2008 Elsevier Ltd. All rights reserved.