function root = secant_demo(fname, x, tol);

% SECANT_DEMO  Secant method (demo)
%
% Input       fname  : Name of function.
%             x      : Vector with two approximations of the root.
%             tol    : Relative error of the found root.
% Output      root   : Position of the root.
%
% Remarks  1) There are functions for which this method fails. So be very 
%             careful using this algorithm!
%          2) Disadvantage of this method is that you need two start values.
%             With a wrong choice of these start values the algorithm will
%             not work properly.

ZoomOn = 0;

if nargin < 3, tol = 1e-4; end
if nargin < 2, x = [0,0.1]; end
if nargin < 1, fname = 'fff'; end


h1 = [];
h2 = [];
clf

lower = min(x);     % can be completely wrong!!
upper = lower+1.2;  % can be completely wrong!!
x = lower:(upper-lower)/999:upper;
y = feval(fname,x);
clf
plot(x,y,x,zeros(size(x)),'--');
axis([min(x),max(x),min(y),max(y)]);
title('Secant Demo (<enter> in command window gives next point)')

x1 = x(1);
x2 = x(2);


while 1,

% Calculate the new approximation based on the secant method

  f1 = feval(fname,x1);
  f2 = feval(fname,x2);
  root = x2 - f2*(x2-x1)/(f2-f1);
  if length(h1),
    delete(h1);
    delete(h2);
  end
  h1 = line([x2,root],[f2,0],'LineStyle',':');
  fnew = feval(fname,root);
  h2 = line([root,root],[0,fnew],'LineStyle',':');
  hold on
  plot(x2,f2,'r.',root,fnew,'r.');

  if ZoomOn,
    minx = min(x2,root);
    maxx = max(x2,root);
    dx = maxx-minx;
    minx = minx-0.1*dx;
    maxx = maxx+0.1*dx;
    miny = min(f2,fnew);
    maxy = max(f2,fnew);
    dy = maxy-miny;
    miny = miny-0.1*dy;
    maxy = maxy+0.1*dy;
    axis([minx,maxx,miny,maxy]);
  end
  pause

  err = abs((root-x2)/root);
  x1 = x2;
  x2 = root;
  if err < tol,
    break;
  end
end

if ZoomOn,
  axis([min(x),max(x),min(y),max(y)]);
end
title('Secant Demo (end)')