In this paper, a higher-order accuracy method is proposed for the solution of time-dependent nature convection problems based on the stream function-vorticity form of Navier-Stokes equations, in which an optimized third-order upwind compact scheme (Opt-UCD3) with high resolution is proposed to approximate the nonlinear convective terms, the fourth-order symmetrical Fade compact scheme is utilized to discretize the viscous terms, the fourth-order compact scheme on the nine-point 2D stencil is used for approximating the stream-function Poisson-type equation and the third-order TVD Runge-Kutta method is employed for the time discretization. To assess numerical capability of the newly proposed algorithm, particularly its spatial behavior, a problem with analytical solution and another one with a steep gradient are numerically solved. Moreover, the nature convection flows in the square cavity with adiabatic horizontal walls and differentially heated vertical walls are also computed for the wide range of Rayleigh numbers (10(3) < Ra < 10(10)). The characteristic parameters such as Nusselt number, velocity, and streamline show excellent agreement with benchmark solutions and some accurate results available in the literature. Additionally, the detailed features of flow phenomena for the higher Rayleigh numbers (10(8) < Ra < 10(10)) are delineated. The results show that the natural convection flow looses stability firstly via a Hopf bifurcation at Ra-c1 to the periodic flow regime, and then undergoes second bifurcation at a critical Rayleigh number Ra-c2 to quasi-periodic flow regime, and eventually transits to turbulent through a further bifurcation. In the periodic regime, there exist at least two branches of solutions. All of the results are agree well with ones in the literature and show the capabilities of the present method to properly simulate the unsteady nature convection problems. (C) 2016 Elsevier Ltd. All rights reserved.