%
%   NLDEMO2_LESSON6.M     Solve nonlinear demo using Newton's method
%
%   This file computes the solution for a 3rd order nonlinear system via
%   Newton's method.  This method requires explicit evaluation of the nonlinear
%   function and the matrix of first partial derivatives (Jacobian matrix).  With 
%   a given guess for x, f(x) and J(x) can be determined, and the solution of a  
%   linear system, J(x) * h = f(x), can be obtained to find the next estimate, 
%   x_i+1 = x_i - h .  This process is continued until convergence (or we hit the 
%   max number of iterations).
%
%   Note that NLDEMO1 and NLDEMO3... solve this same problem using other methods...
%
%   File prepared by J. R. White, UMass-Lowell (last update:  Nov. 2017)
%
      
      clear all;   close all;   nfig = 0;
%
%   set a reasonable initial guess 
      x2 = input(' input a guess for x2: ');  
      if x2 == 3, x2 = 3.001; end   % don't allow x2 = 3 as guess
      x1 = (-7+2*x2*x2)/3;          % satisfies eqn. 1
      x3 = -5/(3-x2)  ;             % satisfies eqn. 3
      xold = [x1 x2 x3]';     
%
%   start iteration loop
      itmax = 30;   it = 0;   tol = 1e-5;   emax = 1;   n = length(xold);
      fprintf(1,'\n Intermediate edit for NLDemo2  \n');
      while emax > tol && it <= itmax
        it = it+1;   
        x = xold;
% 
%   compute function vector using xold
        f = [3*x(1)-2*x(2)*x(2)+7;
              -x(1)*x(2)+3*x(2)-2*x(3)-9;
              -x(2)*x(3)+3*x(3)+5];
%
%   compute Jacobian matrix evaluated at xold
        J = [  3    -4*x(2)     0;
             -x(2)  -x(1)+3    -2;
               0     -x(3)   -x(2)+3];
%
%   compute xnew
        xnew = xold - J\f;
%
%   calc & edit error (intermediate results)
        emax = max(abs((xnew-xold)./xnew));
        fprintf(1,'  it = %3d      max error = %8.3e \n',it,emax);
        fprintf(1,'      xnew         xold     \n');
        for j = 1:n
          fprintf(1,' %10.5f   %10.5f    \n',xnew(j),xold(j));
        end
%
        xold = xnew;  % use current estimate as guess for next iteration
      end
%
%   print final max relative error and iteration count
      fprintf(1,'\n Number of iterations to convergence = %3d\n',it);
      fprintf(1,' Max relative error at convergence =  %8.3e\n',emax);
      if it >= itmax 
        fprintf(1,' ***** WARNING  --  Hit max number of iterations!!!  *****\n');
      end
%
%   end of file