CFD Online Discussion Forums - SIMPLE Algorithm Finite Difference Equations: how to discretize and solve?

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 $(\xi,\eta,t)$ 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

$\frac{\partial u}{\partial t}+\frac{\partial u^2}{\partial\xi}\xi_x+\frac{\partial uv}{\partial\xi}\xi_y-\frac{\partial u}{\partial\xi}\left(\frac{\xi_{xx}+\xi_{yy}}{Re_c} \right )+\frac{\partial u^2}{\partial\eta}\eta_x+\frac{\partial uv}{\partial\eta}\eta_y-\frac{\partial u}{\partial\eta}\left(\frac{\eta_{xx}+\eta_{yy}}{Re_c} \right )$
$=-\frac{\partial p}{\partial\xi}\xi_x-\frac{\partial p}{\partial\eta}\eta_x+\frac{\partial^2u}{\partial\eta\partial\xi}\left[\frac{2\left(\eta_x\xi_x+\eta_y\xi_y\right)}{Re_c}\right]+\frac{\partial^2u}{\partial\xi^2}\left(\frac{\xi^2_x+\xi^2_y}{Re_c}\right)+\frac{\partial^2u}{\partial\eta^2}\left(\frac{\eta^2_x+\eta^2_y}{Re_c}\right)$

Applying the ADI method, I'll first determine the solution for a half time step implicitly in the $\xi$ direction and explicitly in the $\eta$ direction and then a next half time step explicitly in $\xi$ and implicitly in $\eta$ . 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 $(u^{n+1}_{i,j})^2$ and $(uv)^{n+1}_{i+1,j}$ 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: $\frac{\partial u^2}{\partial\xi}$ and $\frac{\partial uv}{\partial\xi}$ . I've seen in several books where they let $F=\frac{u^2}{2}$ 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!