This study aims to implement a more efficient and stable numerical simulation of the hydraulic transient in a complex pumping station system. A finite volume method (FVM) Godunov scheme was established to simulate the simple pipeline and complex pumping station system. The FVM was then introduced to discretize the mathematical models, while the Riemann solver was selected to solve the discrete flux. The MUSCL-Hancock method was utilized to reconstruct the numerical data at the interface of control volumes. The higher numerical accuracy and stability were realized in the Godunov scheme, compared with the frequently-used method of characteristics. Meanwhile, the MINMOD slope limiter was used to avoid false oscillation. The boundary processing of the dual virtual unit was then presented for the second-order accuracy of both the computational region and the boundary, particularly for the simpler computation. The simulation of the improved model was in good agreement with the exact solution and the classical examples. The sensitivity analysis was also performed on the Courant and grid number. Furthermore, a more accurate, stable, and efficient performance was achieved in the second-order Godunov scheme, compared with the method of characteristics. More importantly, there was more outstanding attenuation with the decrease of the Courant number for a simple pipeline system. The computation time of the second-order Godunov scheme was 0.017 s at the same accuracy, compared with the method of characteristics (0.227 s). Consequently, a more stable and efficient performance was achieved in the second-order Godunov scheme. In the actual pumping system with the multiple-characteristics pipe structure, the second-order Godunov scheme required only a slight reduction in the Courant numbers, indicating the simple and convenient way for high numerical accuracy. Once the method of characteristics was used to calculate the hydraulic transition of the pumping station, the Courant number in the pipeline was less than 1 at the same length or wave velocity of the pipeline. By contrast, the Courant number was 0.72-0.76 in this case, indicating a very different simulation from the actual. Therefore, it is necessary to adjust the local pipeline length or wave velocity for the condition that the Courant number was 1. The tedious operation can lead to calculation errors, due to the change in pipeline characteristics. The accuracy can be improved but with less computational efficiency, if the wave velocity remained unchanged to increase the number of computational grids. In the method of characteristics, the number of grids can properly improve the accuracy of the calculation but with the doubled computation time, when the Courant number was less than 1. In the second-order Godunov scheme, there was little effect of grid number on the accuracy of the calculation but with the longer calculation time, whether the Courant number was equal to or less than 1. Therefore, a finer grid was preferred in the method of characteristics for the same accuracy requirements, when the Courant number was less than 1 in the transient process of the simulated pump system. Therefore, the second-order Godunov scheme can accurately simulate the process lines of rotational speed, discharge, and outlet pressure parameters during the hydraulic transient of the pump system. Anyway, the second-order Godunov scheme can be expected to effectively improve the efficiency, stability, and accuracy of hydraulic transient simulation of traditional pumping station systems.