function [x,A,b] = gauss_elim_step_by_step(A,b)
% Solves A*x = b by Gaussian elimination showing every intermediate step.
% b must be a column vector

n = length(b);
A
b

% Elimination phase
for k = 1:n-1                   
    for i= k+1:n
        if A(i,k) ~= 0
            mult = A(i,k)/A(k,k);            
            A(i,k:n) = A(i,k:n) - mult*A(k,k:n); 
            b(i)= b(i) - mult*b(k);
            % elim_mat = [A b [n k i]']
            elim_mat = [A b]
            % A
            % b
         end
     end
end
A
b
b_copy = b;

% Back substitution phase
for k = n:-1:1                  
    b(k) = (b(k) - A(k,k+1:n)*b(k+1:n))/A(k,k);
    % backsub_mat = [A b [n k i]']
    backsub_mat = [A b]
    x = b;
end

b = b_copy;

end
