function [x,dx] = gauss_seidel(A,b,x,eps,N)
%
%  The function gauss_seidel applies the method of Gauss-Seidel
%  to solve A*x = b.
%  On entry:
%    A       coefficient matrix of the linear system;
%    b       right-hand side vector of the linear system;
%    x       initial approximation;
%    eps     accuracy requirement, stop when norm(dx) < eps;
%    N       maximal number of steps allowed.
%  On return:
%    x       approximation for the solution of A*x = b;
%    dx      vector last used to update x, 
%            if success, then norm(dx) < eps.
%
%  Example:
%    Q = orth(rand(5)); D = diag(rand(1,5));
%    A = Q*D*Q'; z = rand(5,1); b = A*z;
%    x = z + rand(5,1)*1e-4;
%    [x,dx] = gauss_seidel(A,b,x,1e-8,50)
%