% Poisson solve
  pold   = p;
  change = 2*eps;
  it     = 1;
  while (change > eps)
      
      pold = p;
      % boundary condition
      %Example p(2:N-1 ,1    ) = p(2:N-1 ,2    )   + 2.0*v(2:N-1 ,3    )/(Re*dx); % bottom 
      
      % NOTE: All boundary condtions must be programmed here (i.e. for all walls!!!!!)
      
      % SOR
      for i=2:N-1
        for j=2:N-1
          p(i,j) = 0.25*(p(i-1,j)+p(i,j-1)+p(i+1,j)+p(i,j+1) - Q(i-1,j-1)*dx^2);
          p(i,j) = pold(i,j) + rf*(p(i,j)-pold(i,j));
        end
      end
    
      
      pmax = max(abs(pold(:)));
      if (pmax == 0) 
          pmax = 1.0;
      end
      change =  max(abs( (pold(:)- p(:))/pmax ));
      it = it + 1;
      if (it > itmax) 
          break;
      end    
  end