%function lu_factorization ( )

  n = 3;

  a = [ 2, -1,  2; ...
        4,  1,  9; ...
        8,  5, 24 ]

  b = [ 0, -5, -16 ]'

%
%  Compute the L U factors.
%
  l = zeros ( n, n );
  u = zeros ( n, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  L | U | L * U\n' );
  fprintf ( 1, '\n' );
  [l,u,l*u]
  pause

  for k = 1 : n

    for j = k : n
      u(k,j) = a(k,j) - l(k,1:k-1) * u(1:k-1,j);
    end

    l(k,k) = 1.0;
    for i = k + 1 : n
      l(i,k) = ( a(i,k) - l(i,1:k-1) * u(1:k-1,k) ) / u(k,k);
    end

    fprintf ( 1, '  L | U | L * U after step %d\n', k );
    [ l, u, l * u ]
    pause

  end


