function roots ()
% test program for root finding routine.
  fprintf('  %8s  %16s  %16s  %6s\n', 'function', 'xend', 'error', 'nmbf');
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func0', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func0', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func1', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func1', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func2', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func2', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func3', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func3', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func4', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func4', xend, xend-sqrt(2), nmbf);
  a0 = 1; b0 = 2;
  [xend, nmbf] = findroot ('func5', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func5', xend, xend-sqrt(2), nmbf);
  a0 = -1; b0 = 2;
  [xend, nmbf] = findroot ('func6', a0, b0);
  fprintf('  %8s  %16.7e  %16.7e  %6d\n', 'func6', xend, xend-0, nmbf);

  function [xend, nmbf] = findroot (fname, a0, b0)
%
% findroot: locate a root or a sign change of a given function in
% a given interval.
% Inputs:
%   f: double function of one double argument.
%   a0, b0: double, bounding bracket for the root.
%   Must have signum(f(a0))*signum(f(b0))<=0.
% Returns: 
%   xend: root of f in the interval defined by (a0,b0).
%   nmbf: number of calls to f made for this call to findroot.
  ...

  function s = signum(x)
% returns the sign of x; returns 0 if x==0.
  if x<0
    s = -1;
  elseif x==0
    s = 0;
  elseif x>0
    s = 1;
  else
    s = 0/0;
  end

  function y = func0(x)
% test function for rootfinding routine.
% func0 is smooth and has a simple root at x=sqrt(2).
  y = (x-sqrt(2))*sqrt(1+x^2);

  function y = func1(x)
% test function for rootfinding routine.
% func1 is smooth and has a triple root at x=sqrt(2).
  y = (x-sqrt(2))^3;

  function y = func2(x)
% test function for rootfinding routine.
% func2 is smooth and has a simple root at x=sqrt(2); the
% asymptotic behavior may pose problems for a Newton-like method.
  y = atan(x-sqrt(2));

  function y = func3(x)
% test function for rootfinding routine.
% func3 is smooth and has a simple root at x=sqrt(2); the
% asymptotic behavior may pose problems for a Newton-like method.
  y = atan(256*(x-sqrt(2)));

  function y = func4(x)
% test function for rootfinding routine.
% func4 is discontinuous with sign change at x=sqrt(2).
  if x<0|x*x<2
    y = -1;
  else
    y = 1;
  end

  function y = func5(x)
% test function for rootfinding routine.
% func5 has a vertical asymptote and sign change at x=sqrt(2).
  t = (x-sqrt(2))*sqrt(1+x^2);
  if t~=0
    y = 1/t;
  else
    y = 0;
  end

  function y = func6(x)
% test function for rootfinding routine.
% func6 is discontinous with sign change at 0.
  if x<0
    y = -1;
  else
    y = 1;
  end