%function gauss_elimination ( )

  n = 3;

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

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

%
%  Forward elimination.
%  For each row I, let E = A(I,I),
%    For each later row I2, let D = A(I2,I),
%      subtract E/D times row I from row I2
%      subtract E/D times B(I) from B(I2).
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A | B\n' );
  fprintf ( 1, '\n' );
  [a,b]
  pause

  for i = 1 : n - 1
    d = a(i,i);
    for i2 = i + 1 : n
      e = a(i2,i);
      a(i2,i:n) = a(i2,i:n) - ( e / d ) * a(i,i:n);
      b(i2) =     b(i2)     - ( e / d ) * b(i);
    end
    fprintf ( 1, '  A | B after elimination step %d\n', i );
    [ a, b ]
    pause
  end
%
%  Back substitution.
%  Solve for variable I.
%  Subtract A(I2,I) * X(I) from preceding equations.
%  Zero out A(I2,I).
%
  x = zeros(n,1);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  A | X | B\n' );
  fprintf ( 1, '\n' );
  [a, x, b]
  pause
  for i = n : - 1 : 1
    x(i) = b(i) / a(i,i);
    for i2 = 1 : i - 1
      b(i2) = b(i2) - a(i2,i) * x(i);
      a(i2,i) = 0.0;
    end
    fprintf ( 1, '\n' );
    fprintf ( 1, '  A | X | B after backsolve step %d\n', i );
    fprintf ( 1, '\n' );
    [ a, x, b ]
    pause
  end


