%  Demo_newton.m
%
%  Solve f(x) = exp(x) - 5 = 0 using Newton's method.

  disp(' Use Newton Method to solve f(x) = exp(x) - 5 = 0.')
  format long;      %  Display 15 digits instead of usual 5.
  x_true = log(5)
  x0 = 3; %%%input(' Initial guess for solution (should be > x_true): ');
  maxiter = 6; %%%input(' Max. no. of Newton iterations = ');
  Ftol = 1e-6; %%%input(' Absolute stopping tolerance, Ftol = ');

%  Save initial guess and initial absolute error in vectors xvec and evec.

  xvec = x0;
  evec = abs(x0 - x_true);

%  Set up vector of t-values between x_true and the initial guess x0. 
%  Evaluate f at these t-values and plot it.  The command "hold on" holds 
%  the plot.  Subsequent plots are written over this one.
  
  t = [x_true:(x0-x_true)/100:x0]';
  F = exp(t) - 5;
  figure(1)
    clf     %  Clear previous graph
    zerovec = zeros(size(t));  %  Dummy variable. Plot it to make legend work.
    plot(t,F,'--', t,zerovec,'-', x_true,0,'*', x0,0,'o')
    grid
    xlabel('x axis')
  hold on

%  Newton initialization.

  k = 1;                                   % Iteration count.
  xk = x0;                                 % Last approx. to solution.
  Fk = exp(xk) - 5;                        % F(xk)
  Fkderiv = exp(xk);                       % F'(xk)
   
%  Newton iteration

  while k < maxiter & abs(Fk) > Ftol,
    
    %  Compute new approximation x_new to solution and plot line segment 
    %  between (xk,Fk) and (x_new,0).
    
    sk = -Fk / Fkderiv;
    x_new = xk + sk; 
    plot(t,F,'--', [xk,x_new],[Fk,0],'-')
    legend('F(x)','Lk(x)','Location','NorthWest')
    plot( x_new,0,'o')  %  Put circle at new solution.
    disp(' Hit any key to continue.'); pause
    
    %  Reset xk, F(xk), and derivative F'(xk).  
    %  Also increment iteration count k.
    
    xk = x_new;
    Fk = exp(xk) - 5;
    Fkderiv = exp(xk);
    k = k + 1;
    
    %  Save approximate solution and error.
    
    xvec(k) = xk;
    evec(k) = abs(xk - x_true);
  end

%  Display vectors of approximates solutions and errors.

  fprintf('\n Iter   Approx       Error\n')
  k_last = k;
  for k = 1:k_last
    fprintf(' %2.0f   %10.8f   %6.4e\n', k, xvec(k), evec(k));
  end

  figure(2)
    semilogy(evec,'o-')
    xlabel('Newton Iteration')
    title('Absolute error of approximation vs Newton iteration')
    grid