function [x,Bind,Nind,sol] = simplex_step(c,A,b,Bind,Nind,stop)
if stop, keyboard; end
cb = c(Bind);
cn = c(Nind);
B = A(:,Bind);
N = A(:,Nind);
sol = 0;

xb = B\b;
xn = 0;

lam = B'\(cb);

sn = cn - N'*lam;

if all(sn >= 0)
    x(Bind) = xb;
    x(Nind) = xn;
    sol=1;
    return
end

% use the Dantzig index
[val,qn] = min(sn);
q = Nind(qn);

assert(all(A(:,q) == N(:,qn)));

d = B\A(:,q);

if all(d <= eps(1))
    sol = -1;
    x = zeros(size(A,2),1);
    return
end

xq = xb./d; xq(d < eps(i)) = Inf;
[xq pb] = min(xq);
p = Bind(pb);

Bind = setdiff(Bind,p); Nind = union(Nind,p); % remove p from B and add to N
Nind = setdiff(Nind,q); Bind = union(Bind,q); % remove q from N and add to B

% Set the next iteration
B = A(:,Bind);
N = A(:,Nind);
xb = B\b;
x(Bind) = xb;
x(Nind) = 0;