I'm attempting to plot velocity and pressure profiles for air flow with

around a NACA airfoil using the SIMPLE procedure to solve the incompressible Navier-Stokes equations. So far I've generated an elliptic grid and transformed the equations pertinent to the SIMPLE procedure to the computational

plane. Where I'm completely confused is how to discretize and solve velocities and pressure. From what I've read, which is a lot of CFD books, I'll be solving the momentum equations for the velocities by the alternating direction implicit (ADI) method. Does this apply to all steps of the SIMPLE procedure (i.e. solving for velocities, solving the pressure correction, and solving velocity corrections)? Or solely the solution to the momentum equations?
The part that really has me confused is how to treat the nonlinear terms. For example, the 2-D x-momentum Navier-Stokes equation is
Applying the ADI method, I'll first determine the solution for a half time step implicitly in the

direction and explicitly in the

direction and then a next half time step explicitly in

and implicitly in

. When I disretize the momentum equations, solving via ADI will leave me with a multidiagonal system. But I'm lost on exactly how to solve because I'll have terms such as

and

as part of the solution vector. I'm actually lost on how to organize the equation to set up as a matrix equation in the first place. I've written out the equations, which has spanned many pages, yet still do not see how to proceed.
My last question is how to discretize the following two terms:

and

. I've seen in several books where they let

but I'm unsure how that gets incorporated into the solution. As for both the nonlinear terms I mentioned, will I basically have

and

values for each point in the domain? Like I said, I'm really confused. I'm also a newbie to CFD so please take it easy on me. I appreciate any advice. Thanks!