It is essential for a Navier–Stokes equations solver based on a projection method to be able to solve the resulting Poisson equation accurately and efficiently. In this paper, we present numerical solutions of the 2D Navier–Stokes equations using the fourth-order generalized harmonic polynomial cell (GHPC) method as the Poisson equation solver. Particular focus is on the local and global accuracy of the GHPC method on non-uniform grids. Our study reveals that the GHPC method enables the use of more stretched grids than the original HPC method. Compared with a second-order central finite difference method (FDM), global accuracy analysis also demonstrates the advantage of applying the GHPC method on stretched non-uniform grids. An immersed-boundary method is used to deal with general geometries involving the fluid–structure interaction problems. The Taylor–Green vortex and flow around a smooth circular cylinder and square are studied for the purpose of verification and validation. Good agreement with reference results in the literature confirms the accuracy and efficiency of the new 2D Navier–Stokes equation solver based on the present immersed-boundary GHPC method utilizing non-uniform grids. The present Navier–Stokes equations solver uses second-order central FDM and Quadratic Upstream Interpolation for Convective Kinematics scheme for the discretization of the diffusion term and advection term, respectively, which may be replaced by other higher-order schemes to further improve the accuracy.