A two dimensional thermal-hydraulic analysis of a natural circulation experiment has been performed to evaluate the effectiveness of a higher order finite difference method for solving the Navier-Stokes and the energy equations. In the method, the convection terms appearing in each equation are solved by the Method of Characteristics using the third order Lagrange type polynomial as the interpolation function, and an iterative procedure is applied to solve the time derivative terms of each equation stably with second order accuracy. The analytical results have been compared with an experiment in which the temperature and the velocity distributions in the plenum region were measured with their fluctuations, and it was shown that the higher order finite difference method could simulate natural convection phenomena involving fluctuations well.