A high-order solver for the blood flow is developed and analyzed using a two-dimensional backward-facing step. In the first part, a Newtonian steady code to solve the incompressible Navier–Stokes (N–S) equations has been developed. The accuracy of the code is verified by comparing the results to the experimental results. An exact projection method/fractional-step scheme is used to solve the incompressible N–S equations. Convective terms of the N–S equations are solved using fifth-order WENO spatial operators, and for the diffusion terms, a sixth- order compact central difference scheme is employed. The third-order Runge–Kutta (R–K) explicit time-integrating scheme with total variation diminishing (TVD) is adopted for time discretization. In the second part, the pulsatile behavior of the Newtonian blood flow has been added to the initial program. Thirdly, the numerical code has been extended to include the steady and pulsatile effects in non-Newtonian blood flow. Finally, a practical example of bend tube has been analyzed by extending the two-dimensional code to 3D and the results are compared to already published data.