This paper is concerned with the development of numerical research for the 2D vorticity-stream function formulation and its application in vortex merging at high Reynolds numbers. A novel numerical method for solving the vorticity-stream function formulation of the Navier-Stokes equations at high Reynolds number is presented. We implement the second-order linear scheme by combining the finite difference method and finite volume method with the help of careful treatment of nonlinear terms and splitting techniques, on a staggered-mesh grid system, which typically consists of two steps: prediction and correction. We show in a rigorous fashion that the scheme is uniquely solvable at each time step. A verification algorithm that has the analytical solution is designed to demonstrate the feasibility and effectiveness of the proposed scheme. Furthermore, the proposed scheme is applied to study the vortex merging problem. Ample numerical experiments are performed to show some essential features of the merging of multiple vortices at high Reynolds numbers. Meanwhile, considering the importance of the inversion for the initial position of the vorticity field, we present an iteration algorithm for the reconstruction of the initial position parameters.