backward_substitution()

function x = backward_substitution(U, b)
  n = length(b);
  x = zeros(n, 1);
  x(n) = b(n)/U(n, n);
  for i = n-1:-1:1
    sum = 0;
    for j = i+1:n
      sum = sum + U(i, j) * x(j);
    endfor
    x(i) = (b(i)-sum)/U(i, i);
  endfor
endfunction