The development of numerical methods for the incompressible Navier-Stokes equations received much attention in the past 50 years. Finite element methods emerged given their robustness and reliability. In our work, we choose the P2-P1 finite element for space approximation which gives 2nd-order accuracy for velocity and 1st-order accuracy for pressure. Our research focuses on the development of several high-order semi-implicit time-stepping methods to compute unsteady flows. The methods investigated include backward difference formulae (SBDF) and defect correction strategy (DC). Using the defect correction strategy, we investigate two variants, the first one being based on high-order artificial compressibility and bootstrapping strategy proposed by Guermond and Minev (GM) and the other being a combination of GM methods with sequential regularization method (GM-SRM). Both GM and GM-SRM methods avoid solving saddle point problems as for SBDF and DC methods. This approach reduces the complexity of the linear systems at the expense that many smaller linear systems need to be solved. Next, we proposed several numerical improvements in terms of better approximations of the nonlinear advection term and high-order initialization for all methods. To further minimize the complexity of the resulting linear systems, we developed several new variants of grad-div splitting algorithms besides the one studied by Guermond and Minev. Splitting algorithm allows us to handle larger flow problems. We showed that our new methods are capable of reproducing flow characteristics (e.g., lift and drag parameters and Strouhal numbers) published in the literature for 2D lid-driven cavity and 2D flow around the cylinder. SBDF methods with grad-div stabilization terms are found to be very stable, accurate and efficient when computing flows with high Reynolds numbers. Lastly, we showcased the robustness of our methods to carry 3D computations.