gauss

function [ x, U ] = gauss( A, b )
%Gaussian Elimination algorithm mimicking
%algorithm 2.3 in the textbook.

n = length(b); %how many row/columns our matrix is
m = zeros(n, 1); %placeholders
x = zeros(n,1); %answer matrix

for k = 1 : (n - 1)
    if (A(k,k) == 0)
        error('Pivot is zero: breaking');
    end

    %Compute a column of M
    for i = (k + 1) : n
        m(i) = (A(i,k) / A(k,k));
    end

    %Compute A and b after pivot operations
    for i = (k + 1) : n
       for j = (k + 1) : n
           A(i,j) = A(i,j) - m(i) * A(k,j);
       end
    end

    for i = (k + 1) : n
       b(i) = b(i) - b(k) * m(i);
    end
end

U= triu(A); %uper triangle of A

%Backward substitution
x(n) = b(n) / A(n,n);
for k = n-1 : -1 : 1;
    b(1:k) = b(1:k)  -x(k+1) * U(1:k,k+1);
    x(k) = b(k) / U(k,k);
end;
end