gauss

function [ A,b,xn] = gauss( A,b )
n = size(A,1);

for col=1:n-1
    for i=1+col:n
        if A(i, col) == 0
            continue
        end
        mul = A(i,col)/A(col,col);
        A(i,:) = A(i,:) - A(col,:)*mul;
        b(i,:) = b(i,:) - b(col,:)*mul;
    end
end

xn = zeros(n,1);
xn(n) = b(n,1)/A(n,n);
for i=n-1:-1:1
    sum = 0;
    for xi=1:n
        sum = sum + xn(xi)*A(i,xi);
    end
    if sum == 0
        xn(i) = 0;
    else
        xn(i) = (b(i)-sum)/A(i,i);
    end
end
end