display('Number of equations:');
n = 3

display('Coefficient matrix:');
A = magic(n)

display('Inhomogenity:');
b = [1; 2; 3]

% Create zeros below the diagonal for each column j
for j = 1:n
    % Create zeros below the diagonal in column j for each row i
    for i = j+1:n
        % Compute the factor with which we have to multiply the j-th
        % row when adding it to the i-th row.
        c = A(i, j) / A(j, j);
        % Just for educational purposes. This particular A(i, j)
        % becomes zero and will never be read again in the remainder
        % of this algorithm.
        A(i, j) = 0;
        % Subtract c * row j from row i.
        for k = j+1:n
            A(i, k) = A(i, k) - c * A(j, k);
        end;
        % Also perform the same manipulation on the inhomogenity.
        b(i) = b(i) - c * b(j);
    end
end

% Backward substitution. Begin with computing the solution
% in the last row and work backwards until the first row.
for i = n:-1:1
    % Solve for variable x(i). Try compute this manually to see
    % what's happening here.
    for j = i+1:n
        b(i) = b(i) - A(i, j) * b(j);
    end
    % Divide by coefficient for variable x(i) in row i.
    b(i) = b(i) / A(i, i);
end

display('Right upper diagonal matrix:');
A

display('Solution:');
b
