%Gauss-Jordan algorithm with pivoting
function [x,AB]=gaussjordan(A,B)
NA=size(A,2);
[NB1,NB]=size(B);
if NB1~=NA
    error('A and B must have compatible dimensions')
end
AB=[A B];
E=eye(NA);
for j=1:NA
    for i=1:NA
        if AB(j,j)~=0
            E(i,j)=-AB(i,j)/AB(j,j);
            AB=E*AB;
            E=eye(NA)
        else
            [maxx,k]=max(abs(AB([j:NA],j)));
            AB([j k+(j-1)],:)=AB([k+(j-1) j],:);
            E(i,j)=-AB(i,j)/AB(j,j);
            AB=E*AB;
            E=eye(NA)
        end
    end
end
E=eye(NA);
for i=1:NA 
    E(i,i)=1/AB(i,i);
end
AB=E*AB;
x=AB(:,end)