
clc,clear

AA = [0.001 2 3;-1 3.712 4.623;-2 1.072 5.643];
b = [1;2;3];

A = [AA b]; % 增广矩阵

for n = 1:size(A,1)
    
    %% 寻找列主元
    m = find( abs(A(n:end,n)) == max(abs(A(n:end,n))) );

    if A(m+n-1,n)==0
        
        disp('列主元有0元素，Gauss列主元消去法失败！')
        break;
        
    elseif m ~= 1
        
        B = A(m+n-1,n:end);
        A(m+n-1,n:end) = A(n,n:end);
        A(n,n:end) = B;
        
    end
    
    M = A(n,n:end)/A(n,n); % 计算 m_ik
    
    %% Gauss 消去（使用Kronock积Kron(A,b)，可去掉循环）
    A(n+1:end,n:end) = A(n+1:end,n:end) - kron(A(n+1:end,n),M); 
    
end

%% 求解
x = zeros(size(A,1),1);

for n = size(A,1):-1:1
    
    if n == size(A,1)
        
       x(n) = A(n,end)/A(n,n);
        
    end
   
    x(n) = (A(n,end)-A(n,n+1:end-1)*x(n+1:end))/A(n,n);
    
end 

fprintf('方程组的解为：[%g, %g, %g]\n\n',x(1),x(2),x(3));
