function [r, xk] = secant(f, x0, x1, tol, max_iter)
% SECANT find the zero of function f by sceant method.
%
% [r, xk] = secant(f, x0, x1, tol, max_iter)
%
% f        - the function
% x0, x1   - initial approximations
% tol      - tolerance
% max_iter - maximum number of iterations
%
% r        - the approximation of the root
% xk       - successive approximations of the root

if nargin < 4
    tol = 1e-6;
end
if nargin < 5
    max_iter = 100;
end

xk    = zeros(1, max_iter+1); 
xk(1) = x0;
xk(2) = x1;

for k=2:max_iter
    % x2 = x1 - f(x1)(x1-x0) / (f(x1) - f(x0));
    d = f(xk(k)) - f(xk(k-1));
    if abs( d ) < tol
        break;
    end

    xk(k+1) = ( xk(k-1) * f(xk(k)) - xk(k) * f(xk(k-1)) ) / d;

    if abs(f(xk(k+1)) - f(xk(k)) ) < tol
        break;
    end
end

if k >= max_iter
    warning('Max. no. of iterations exhausted without reaching required tolerance.');
end

r = xk(k);
xk = xk(1:k);    % trim to the actual number of iterations performed
return

