Natural circulation refers to the phenomenon of fluid circulation driven by the driving force formed by the density difference and height difference. Natural circulation is very common in nuclear power system, which has an important impact on reactor safety. Therefore, accurate simulation of natural circulation phenomenon is one of the most important research contents in current reactor thermal hydraulic analysis. The low order difference method is widely employed to discretize basic conservation equations for current main thermal-hydraulics analysis code such as RELAP5, CATHARE, TRACE, etc. However, the low order difference scheme has the characteristics of large numerical truncation error and low calculation accuracy when simulating natural circulation problem. In this paper, the high-order accuracy numerical algorithm for two-fluid two-pressure model is used to simulate Welander natural circulation problem. Numerical results imply that the higher order difference scheme can accurately predict the unstable boundary and chaotic behavior with fewer grids, while the lower order difference scheme may get wrong results when the mesh is not fine enough, that is, the unstable mode of natural circulation can be predicted as stable mode. Therefore, high-order accuracy numerical algorithm could prevent excessive numerical diffusion effectively and improve the prediction accuracy, which demonstrates the advantage of using high-order schemes.