Thanks again for the help! so you are saying that after solving the momentum equations yielding u* and v*, that i should apply the outlet boundary condition to these guessed velocity fields u* and v*? That makes sense, as that velocity cv (at i=NI), is the location of the east face of the pressure control volume at I=NI-1, i.e. the last pressure cv in the domain.

On one quick note about the initial guesses, which initialise the solution... at the moment, I have a guessed pressure field p* which is a linear distribution from west to east, i.e. 10-8-6-4-2-0, and I have also set the intial velocity field u* equal to the inlet vlaue, and v*=0. Anderson suggests putting in a velocity spike in a v control volume to create a 2D flow also.
If i set the u velocity field to zero (i=3 onward to i=NI-1), my outlet b.c returns NaN, and the solution fails due to this piece of code:

u(nX,2:nYp-1)=u(nX-1,2:nYp-1); % First extrapolate values from the domain

u(nX,2:nYp-1)=u(nX-1,2:nYp-1)*sum(u(2,2:nYp-1))/sum(u(nX,2:nYp-1)); % This multiplier is the mass in divided by the mass out, which ensures continuity

So I guess the question is, do i assign velocity values that correspond with p*, or is it just sufficient that the u-velocities be non-zero?

This is all much appreciated and I am making more progress than I have in a long time! Best Regards, Michael

Yeah you should ensure continuity on outlet every time before solving pressure correction.
As for NaN problem, I guess, because you applied 0 to all u velocity as initial guess, so when you apply
actually you are dividing by zero,or sth very near zero.
As far as I know...pressure initial guess is not so crucial problem.
good luck

I assume then that this is what you mean:

%% Apply Boundary Conditions
boundary_conditions();
%% Solve for u*
[ustar, diJ]=u_momentum();

%% Check continuity again
% u(nX,2:nYp-1)=ustar(nX-1,2:nYp-1); % First extrapolate values from the domain
% u(nX,2:nYp-1)=ustar(nX-1,2:nYp-1)*sum(u(2,2:nYp-1))/sum(u(nX,2:nYp-1));
%% Solve for v*
[vstar, dIj]=v_momentum();

%% Solve for p'
[pdash] = pressure_correction(diJ, dIj);


Where now I have used the values of u* along i=NI-1, whereas in the boundary conditions, I had used the values of the initial guess.

Thank you very much for your notes!
My cavity code wasn't working till I used your suggested relaxation factors! And now it is working like a charm!!