% Programming and Numerical Methods
% MSc in Computational and Applied Physics

% Code 17: Secant method for solving f(x) = 0 in [a,b]
% A. Meseguer (2007)

% On input
% 1. [a,b]: interval (it asumes that f(a)f(b) < 0)
% 2. tol:   tolerance so that abs(bk-ak)/2 < tol
% 3. itmax: stop criteria for exceeding number of iterations
% 4. fun:   function to be solved

% On output
% 1. xk:   resulting iteration points
% 2. resd: resulting residuals iteration
% 3. it:   number of required iterations to satisfy tolerance 

  function [xk,resd,it] = secant(a,b,tol,itmax,rootfun)

% We compute two first points wth the chord method
  fa = feval(rootfun,a) ; fb = feval(rootfun,b) ; q = (fb-fa)/(b-a);
  x0 = b ; x1 = x0 - fb/q ;
  f0 = fb ; f1 = feval(rootfun,x1) ;

  xk = [x0 x1] ; resd = [abs(f0) abs(f1)] ; it = 2 ; 
  tolk = .5*abs(b-a) ;

  while it < itmax & tolk > tol
   q = (f1-f0)/(x1-x0) ; 
   xk = [xk , xk(end) - f1/q] ;  
   x0 = xk(end-1) ; x1 = xk(end) ;
   f0 = f1 ; f1 = feval(rootfun,x1) ;
   resd = [resd abs(f1)] ; tolk = 0.5*abs(xk(end)-xk(end-1)) ; 
   it = it + 1 ;
  end





