%% Gaussian Elimination algorithm with pivoting
function [x,AB]=gauss2(A,B)
NA=size(A,2);
[NB1,NB]=size(B);
if NB1~=NA
    error('A and B must have comppatible dimensions')
end
AB=[A B];
E=eye(NA,NA);
for j=1:NA-1
    for i=j+1:NA
        if AB(j,j)~=0
            E(i,j)=-AB(i,j)/AB(j,j);
            AB=E*AB;
            E=eye(NA,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,NA);
        end
    end
end
% backward substitutation for a upper-triangular matrix equation
x=zeros(1,NA);
sum=0;
for i=NA:-1:1
    for j=i+1:NA
        sum=sum+x(1,j)*AB(i,j);
    end
    x(1,i)=(AB(i,NA+1)-sum)/AB(i,i);
    sum=0;
end
