In the nemurical fluid analysis, high accuracy and upwinding are necessary for the convection term of the Navier-stokes equation because this term is non-linear. In general, bilinear element is used in FEM (Finite Element Method) and the discretization by this method corresponds to the 2nd-order center difference. As a result, we can't obtain adequate accuracy for high Reynolds number flow. Although many upwind methods have been proposed, we use CIP (Cubic Interpolated Pseudoparticle)-FEM in the present work, and analyze natural convection in a square cavity. This problem has been analyzed by many investigators. However, it is difficulut to analyze this flow at high Rayleigh numbers Ra≥10⁹ because numerical divergences and vibrations occur in energy equation. We analyze natural convection from 10³≤Ra≤10⁹ and show that the present scheme analyze this flow stably and is effective by comparing with other numerical results.