Implementing Simplex Method infinite loop

I am trying to implement a simplex algorithm following the rules I was given at my optimization course. The problem is

min c'*x    s.t.
    Ax = b
    x >= 0

All vectors are assumes to be columns, ' denotes the transpose. The algorithm should also return the solution to dual LP. The rules to follow are:

Here, A_J denotes columns from A with indices in J and x_J, x_K denotes elements of vector x with indices in J or K respectively. Vector a_s is column s of matrix A.

Now I do not understand how this algorithm takes care of condition x >= 0, but I decided to give it a try and follow it step by step. I used Matlab for this and got the following code.

X = zeros(n, 1);
Y = zeros(m, 1);

% i. Choose starting basis J and K = {1,2,...,n} \ J
J = [4 5 6]  % for our problem
K = setdiff(1:n, J)

% this while is for goto
while 1
    % ii. Solve system A_J*\bar{x}_J = b.
    xbar = A(:,J) \ b

    % iii. Calculate value of criterion function with respect to current x_J.
    fval = c(J)' * xbar

    % iv. Calculate dual solution y from A_J^T*y = c_J.
    y = A(:,J)' \ c(J)

    % v. Calculate \bar{c}^T = c_K^T - u^T A_K. If \bar{c}^T >= 0, we have
    % found the optimal solution. If not, select the smallest s \in K, such
    % that c_s < 0. Variable x_s enters basis.
    cbar = c(K)' - c(J)' * inv(A(:,J)) * A(:,K)
    cbar = cbar'

    tmp = findnegative(cbar)
    if tmp == -1  % we have found the optimal solution since cbar >= 0
        X(J) = xbar;
        Y = y;
        FVAL = fval;
        return
    end

    s = findnegative(c, K)  %x_s enters basis

    % vi. Solve system A_J*\bar{a} = a_s. If \bar{a} <= 0, then the problem is
    % unbounded.
    abar = A(:,J) \ A(:,s)

    if findpositive(abar) == -1  % we failed to find positive number
        disp('The problem is unbounded.')
        return;
    end

    % vii. Calculate v = \bar{x}_J / \bar{a} and find the smallest rho \in J,
    % such that v_rho > 0. Variable x_rho exits basis.
    v = xbar ./ abar
    rho = J(findpositive(v))

    % viii. Update J and K and goto ii.
    J = setdiff(J, rho)
    J = union(J, s)
    K = setdiff(K, s)
    K = union(K, rho)
end

Functions findpositive(x) and findnegative(x, S) return the first index of positive or negative value in x. S is the set of indices, over which we look at. If S is omitted, whole vector is checked. Semicolons are omitted for debugging purposes.

The problem I tested this code on is

c = [-3 -1 -3 zeros(1,3)];
A = [2 1 1; 1 2 3; 2 2 1];
A = [A eye(3)];
b = [2; 5; 6];

The reason for zeros(1,3) and eye(3) is that the problem is inequalities and we need slack variables. I have set starting basis to [4 5 6] because the notes say that starting basis should be set to slack variables.

Now, what happens during execution is that on first run of while, variable with index 1 enters basis (in Matlab, indices go from 1 on) and 4 exits it and that is reasonable. On the second run, 2 enters the basis (since it is the smallest index such that c(idx) < 0 and 1 leaves it. But now on the next iteration, 1 enters basis again and I understand why it enters, because it is the smallest index, such that c(idx) < 0. But here the looping starts. I assume that should not have happened, but following the rules I cannot see how to prevent this.

I guess that there has to be something wrong with my interpretation of the notes but I just cannot see where I am wrong. I also remember that when we solved LP on the paper, we were updating our subjective function on each go, since when a variable entered basis, we removed it from the subjective function and expressed that variable in subj. function with the expression from one of the equalities, but I assume that is different algorithm.

Any remarks or help will be highly appreciated.