function root = newton_demo(fname, x0, tol);

% NEWTON_DEMO  Newton-Raphson method (demo)
%
% Input       fname  : Name of function.
%             x0     : First apoproximation 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 algoritm!
%          2) A disadvantage of this method is that you need the first
%             derivative of the function.

ZoomOn = 0;

if nargin < 3, tol = 1e-4; end
if nargin < 2, x0 = 0; end
if nargin < 1, fname = 'fff'; end

oldroot = x0;
h1 = [];
h2 = [];
clf

lower = x0;      % can be completely wrong!!
upper = x0+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('Newton-Raphson Demo (<enter> in command window gives next point)')

while 1,

% Calculate the new approximation based on the Newton-Raphson method

  [f,df] = feval(fname,oldroot);
  root = oldroot - f/df;
  if length(h1),
    delete(h1);
    delete(h2);
  end
  h1 = line([oldroot,root],[f,0],'LineStyle',':');
  fnew = feval(fname,root);
  h2 = line([root,root],[0,fnew],'LineStyle',':');
  hold on
  plot(oldroot,f,'r.',root,fnew,'r.');

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

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

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