function x = Gauss_Jacobi(FunFcn,x0,nmax,tol,trace)

% Gauss Jacobi method for solving a system of nonlinear equations:

%

% USAGE: X = gauss_ja( 'f_name', x0, max_no, Tol, Trace)  

%            max_no is the maximum number of iterations allowed

%            Tol    is the stopping tolerance

%            Trace  is the option (=1) for printing all iterates.

% For example, X = Gauss_Ja('xyz',[1 1 1]') 

%            solves the following system                          

%          /  sin(x) + y^2 + log(z) - 7 = 0,

%         ||        3*x + 2^y - z^3 + 1 = 0,

%          \              x + y + z - 5 = 0.

%       Here the M-file, xyz.m, contains the following:

%       function q = xyz(p)

%         x = p(1); y = p(2); z = p(3);

%         q(2) = sqrt(7-sin(x)- log(z));

%         q(1) = - (2^y - z^3 + 1)/3;    %%  See also  

%         q(3) = 5 - x - y;           %%%%%  full_new,  nesolve

if nargin<2, help gauss_ja, return, end

if nargin<5 

    trace=0;

end;



if nargin<4

    tol=0.001;

end;



if nargin<3 

    nmax=20;

end;



% these default settings of nmax and tol are arbitrary



x=[];g=[];

x1=x0;



for n=1:nmax

  [g]=feval(FunFcn,x1);

  if norm(g)==Inf, x=Inf; return; end %% Stop failures

  if norm(g)==NaN, x=Inf; return; end %% Stop failures

  if norm(x1-g)<tol

     n=nmax;

  elseif trace==0

        x1=g;x=x1';

     else

        x1=g;x=[x; x1'];

     end

  end

end

