I need to write solution for 2 phase Simplex method. I made solution, but I’m not sure that it works correctly and sometimes answers are not correct.

For example, my code gives me solution:
0 0 0.666667 3.250000
optimal value = 1.416667

Other code gives me diffetent solution:
optimal solution:
0 0 3.5000 0 0 6.5000
objective value: 24.5

If you have any suggestions how to improve it, please, help. Thanks on advance.

C = [3;4;7;-1];
A = [1 3 2 2;
     2 4 1 3];
b = [7; 10];
dn = zeros(4,1);
dv = inf(4,1);

% initialize Jb 
Jb = [3 4];

inv(A(:, Jb))*b;

% dimentions of task 
[m,n] = size(A);
itermax = 10;

% Basis plan
% 1st phase

X = [0; 0; 2/3; 13/4];

for iter = 1:itermax
    % potentials and evaluations
    Jn = setdiff(1:n, Jb);
    Cb = C(Jb);
    u = Cb' * inv(A(:, Jb));
    del = C' - u*A;

    % check optimality criterion
    J0 = [];
    for j = Jn
        if del(j) <= 0 && X(j) == dn(j)
            fprintf('All rightn');
        elseif del(j) >= 0 && X(j) == dv(j)
            fprintf('All rightn');
        else
            J0 = [J0 j];
        end
    end

    if isempty(J0)
        [~, ind] = max(abs(del(J0)));
        J0 = J0(ind);
        
        % build direction l
        l = zeros(n, 1);
        l(J0) = sign(del(J0));
        lb = -inv(A(:, Jb)) * A(:, J0) * sign(del(J0));

        %steps
        theta0 = dv(J0) - dn(J0);
        theta = [];
        for J = Jb
            if l(J) > 0
                theta = [theta (dv(J) - X(J)) / A(:, J)'*lb];
            elseif l(J) < 0
                theta = [theta (dn(J) - X(J)) / l(J)];
            else
                theta = [theta inf];
            end
        end

        [theta_min, ind2] = min([theta theta0]);
        if ind2 <= m
            Jz = Jb(ind2);
        else
            Jz = J0;
        end

        Jb = [setdiff(Jb, Jz) J0];
        X = X + theta_min*lb;
    end
end
fprintf('optimal plan = %fn', X);
fprintf('optimal value = %fn', C'*X);