%% Demonstrate Jacobi iteration on a 2x2 matrix.
%
a = [ 2, -1; ...
     -1,  3 ];

x = [ 2; ...
      5 ];

b = a * x;

x = rand ( 2, 1 );

r = norm ( a * x - b );

fprintf ( 1, '  %g  %g  %g\n', r, x(1), x(2) );
pause

while ( 1 )

  xnew(1,1) = ( b(1) - a(1,2) * x(2) ) / a(1,1);
  xnew(2,1) = ( b(2) - a(2,1) * x(1) ) / a(2,2);

  r = norm ( a * xnew - b );

  x = xnew;

  fprintf ( 1, '  %g  %g  %g\n', r, x(1), x(2) );
  pause

end
