% Code 2: Newton's method for solving f(x) = 0 
% Input:   1. a: initial guess
%          2. tol: tolerance so that abs(x_k+1 - x_k) < tol
%          3. itmax: maximum number of iterations allowed
%          4. fun: function's name 
%          5. dfun: derivative's name 
%          (If dfun = '0', then f'(x) is approximated) 
% Output:  1. xk: resulting sequence 
%          2. res: resulting residuals 
%          3. it: number of required iterations  

function [xk,res,it] = newton(a,tol,itmax,fun,dfun)
  xk = [a] ; fk = feval(fun,xk) ; res = abs(fk) ; it = 0 ; 
  tolk = res(1) ; dx = 1e-8 ;
  while it < itmax & tolk > tol
      if dfun == '0' 
          dfk = (feval(fun,xk(end)+dx)-fk)/dx;
      else
          dfk = feval(dfun,xk(end));
      end
      xk = [xk , xk(end) - fk/dfk] ; fk = feval(fun,xk(end)) ; 
      res = [res abs(fk)] ; tolk = abs(xk(end)-xk(end-1)) ; 
      it = it + 1;
  end
