A high-order immersed boundary method (IBM) for the computation of unsteady, incompressible fluid flows on two-dimensional, complex domains is proposed, analyzed, developed and validated. In the IBM, the equations of interest are discretized on a fixed Cartesian grid. As a result, domain boundaries do not always conform to the (rectangular) computational domain boundaries. This gives rise to 'immersed boundaries', i.e., boundaries immersed inside the computational domain. A new IBM is proposed to remedy problems in an older existing IBM that had originally been selected for use in numerical flow control investigations. In particular, the older method suffered from considerably reduced accuracy near the immersed boundary surface where sharp jumps in the solution, i.e., jump discontinuities in the function and/or its derivatives, were smeared out over several grid points. To avoid this behavior, a sharp interface method, originally developed by LeVeque & Li (1994) and Wiegmann & Bube (2000) in the context of elliptic PDEs, is introduced where the numerical scheme takes such discontinuities into consideration in its design. By comparing computed solutions to jump-singular PDEs having known analytical solutions, the new IBM is shown to maintain the formal fourth-order accuracy, in both time and space, of the underlying finite-difference scheme. Further validation of the new IBM code was accomplished through its application to several two-dimensional flows, including flow past a circular cylinder, and T-S waves in a flat plate boundary layer. Comparison of results from the new IBM with results available in the literature found good agreement in all cases