Research on using phase change material (PCM) suspension to improve the heat transfer and energy storage capabilities of thermal systems is booming; however, there are limited studies on the application of PCM suspension in transient natural convection. In this paper, the implicit finite difference method was used to numerically investigate the transient and steady-state natural convection heat transfer in a square enclosure containing a PCM suspension. The following parameters were included in the simulation: aspect ratio of the physical model = 1, ratio of the buoyancies caused by temperature and concentration gradients = 1, Raleigh number (RaT) = 103−105, Stefan number (Ste) = 0.005−0.1, subcooling factor (Sb) = 0−1.0, and initial mass fraction (or concentration) of PCM particles (ci) = 0−0.1. The results showed that the use of a PCM suspension can effectively enhance heat transfer by natural convection. For example, when RaT = 103, Ste = 0.01, ci = 0.1, and Sb = 1, the steady-state natural convection heat transfer rate inside the square enclosure can be improved by 70% compared with that of pure water. With increasing Sb, the Nusselt number can change nonlinearly, resulting in a local optimal value.