function [x,A] = gauss_elim_outer_piv(A,b) % A should be square % size(A) = [length(x) length(x)] n = size(A,1); p = 1:n; for k = 1:(n-1) rows = k:n; [y,v] = max(abs(A(rows,k))); i = rows(v); p([k i]) = p([i k]); mults = A(p((k+1):n),k) / A(p(k),k); A(p((k+1):n),k:n) = A(p((k+1):n),k:n) - mults * A(p(k),k:n); b(p((k+1):n)) = b(p((k+1):n)) - mults * b(p(k)); end x = zeros(n,1); for k = n:-1:1 sum = b(p(k)) - A(p(k),(k+1):n) * x((k+1):n); x(k) = sum / A(p(k),k); end