function [x,A] = gauss_elim(A,b)
% by Shelvean Kapita
% naive Gaussian elimination
% A is an n by n matrix
% b is an n vector
% solves Ax = b
n = size(A,1);
for j=1:n-1
    if abs(A(j,j)) < eps
        error('zero pivot encountered');
    end
    for i = j+1:n
        mult = A(i,j)/A(j,j);
        for k = j:n
            A(i,k) = A(i,k)-mult*A(j,k);
        end
        b(i) = b(i)-mult*b(j);
    end
end
% back substitution
x = zeros(n,1);
for i = n:-1:1
    for j=i+1:n
        b(i)=b(i)-A(i,j)*x(j);
    end
    x(i)=b(i)/A(i,i);
end