function rt = tan1xRoots(xmin,xmax,nxSearch,tol)
% tan1xRoots  Find the roots of f(x) = tan(x) - 1/x = 0
%
% Synopsis:  rt = tan1xRoots
%            rt = tan1xRoots(xmin)
%            rt = tan1xRoots(xmin,xmax)
%            rt = tan1xRoots(xmin,xmax,nxSearch)
%            rt = tan1xRoots(xmin,xmax,nxSearch,tol)
% 
% Input:  xmin,xmax = range of x over with the roots are to be found
%                     Default: xmin = -4*pi, xmax = 4*pi
%         nxSearch = number of search intervals used by the bracketing
%                    function, brackPlot.  Default: nxSearch = 75
%         tol = relative tolerance on x and f(x) used to terminate the
%               bisection iterations.  Default: tol = 5.0e-12
%
% Output:  rt = vector of roots found on the interval, xmin <= x <= xmax

if nargin<1, xmin=-4*pi;  end
if nargin<2, xmax=4*pi;   end
if nargin<3, nxSearch=75; end
if nargin<4, tol=5.0e-12; end

% -- Use brackPlot to identify intervals that may contain roots
xb = brackPlot('tan1x',xmin,xmax,nxSearch);
ymax = 3;
ylim([-ymax, ymax])  %  Set y-axis limits to clip singularities
nb = size(xb,1);
fprintf('%d brackets found on %6.2f <= x <= %6.2f\n',nb,xmin,xmax)

% -- Loop over the brackets.  Store each root temporarily in r
%    Copy the true roots to the rt vector.  kr is the index into rt
kr = 0;
fprintf('\n   i     r           f(r)\n')
for i=1:nb

  % -- Refine the bracket to isolate the root or singularity
  r = bisect('tan1x',xb(i,:),tol,tol);

  % -- Print the suspected root, and the value of f(x) there
  fr = tan1x(r);
  fprintf('%4d  %8.4f  %12.3e',i,r,fr);
  
  % -- Test the magnitude of f(r) to identify true roots
  %    Assume that abs(fr)<5e-9 for a true root
  if abs(fr)<5.0e-9
    kr = kr + 1;    %  True root has been found, increment counter
    rt(kr) = r;     %  Store the true root
    fprintf('\n')   %  End the printout line with no further message
  else
    fprintf('  Singularity!\n')   %  print a message and end the line
  end
end

end