function x = gauss_elim_outer(A,b)

% A should be square
% size(A) = [length(x) length(x)]

n = size(A,1);

for k = 1:(n-1)
    mults = A((k+1):n,k) / A(k,k);
    A((k+1):n,k:n) = A((k+1):n,k:n) - mults * A(k,k:n);
    b((k+1):n) = b((k+1):n) - mults * b(k);
end

x = zeros(n,1);
for k = n:-1:1
    sum = b(k) - A(k,(k+1):n) * x((k+1):n);
    x(k) = sum / A(k,k);
end