function [ x ] = gauss_elim( A, b )
%%gauss_elim uses Gauss elimination with partial pivoting to solve Ax = b
%%Inputs:  A - A matrix
 %         b - b vector
%%Outputs: x - solution to Ax = b
%%Author: Gavin St. John, University of Florida

%append
U = horzcat(A,b);

%iterate through each column which needs to be 0'd -> create upper
%triangular matrix
for ii = 1:size(U,1)-1
    % figure out which
    [~,m] = max(abs(U(ii:end,ii)));
    m = m+(ii-1);
    U([ii m],:) = U([m ii],:);
    % construct a matrix to subtract from U (UppterTriangularModifier)

    vec = repelem(U(:,ii),1,length(U)); % creates the basis for the matrix
    axx = U(ii,ii); % current pivot value
    R = repmat(U(ii,:),size(U,1),1); % matrix of rows of current column

    utmod = (vec./axx).*R; % puts it all together
    utmod(1:ii,:) = 0; % makes all rows above current pivot entry 0 so we don't affect stuff we've already operated on
    if ii > 1, utmod(:,ii-1) = 0; end % after the first itreration make columns up to the current pivot entry 0
    U = U - utmod; % best part
end

bf = U(:,end); % split appended matrix back up
Af = U(:,1:end-1);
x = Af\bf; % solve backward substitution

end