%----------------------------------------------------------------------
% Usage: [Z, err, fZ] = findzeros(x_start, x_step, x_end, fun)
% Find roots of f(x) = 0 in the interval [x_start, x_end].
%----------------------------------------------------------------------
%   fun:  Function handle (e.g. @sin) or the function's name as string.
% Returns:
%    Z :   The zeros.
%  err :   Values abs(a-b), where a, b are the closest bounds such that
%              f(a) <= 0 <= f(b) or
%              f(a) >= 0 >= f(b) and abs(a-b) <= eps.
%   fZ :   Values f(Z)
%-----------------------------------------------------------------------
% The function f(x) is sampled at equally spaced x's of distance x_step,
% then bisection is used to find the zeros of f(x).
% In order to find _all_ zeros, choose x_step small enough, so that:
%  1. There is at most one zero between two samples.
%  2. For every zero x_0 there are samples f(a), f(b) with
%        a < x_0 < b and
%        sign(f(a)) != sign(f(b)).
% Otherwise only some zeros will be found. If all samples are > 0 or < 0
% respectively, the procedure will fail. Remark: Zeros, that are local
% minima/maxima might be "overlooked". This procedure is especially
% _not_ suitable when all zeros are local minima/maxima.
%-----------------------------------------------------------------------
% Example 1: Finding zeros of sin(x) in [0,10].
%    Z = findzeros(0, 0.5*pi(), 10, "sin")
%-----------------------------------------------------------------------
% Example 2: Finding zeros of bessel functions of the first kind. (The
% procedure is suitable for functions with a single argument only, so
% we need helper functions.)
%    function y = J0(x); y = besselj(0,x); endfunction;
%    function y = J1(x); y = besselj(1,x); endfunction;
%    function y = J2(x); y = besselj(2,x); endfunction;
%    [Z0,er0] = findzeros(0, 0.25*pi(), 100, "J0");
%    [Z1,er1] = findzeros(0, 0.25*pi(), 100, "J1");
%    [Z2,er2] = findzeros(0, 0.25*pi(), 100, "J2");
%-----------------------------------------------------------------------
% Warnings: This might also find "zeros" that in reality are poles:
%   > function y = f(x); y = 1 ./ x; end
%   >[z,r,fZ] = findzeros(-1,0.3,1,"f")
%   warning: division by zero
%   z = -4.9407e-324
%   r =  4.9407e-324
%   fZ = -Inf
% Also, precisely "hitting" the zeros with the samples will lead to
% strange results / overlooking of zeros. The following zeros of sin(x)
% are only found due to roundoff errors.
%   > [z,r] = findzeros(0,2*pi(),100,"sin")
%   z =
%       0.00000
%      81.68141
%      84.82300
%   r =
%      4.9407e-324
%      1.4211e-014
%      1.4211e-014
%-----------------------------------------------------------------------
% Tested with Octave 3.2.4; public domain; no warranties whatsoever; not
% built for speed; don't put your life on it; don't sue me :)
% Version: 2010-04-14, S. Reiser (s.reiser@tu-bs.de).
%-----------------------------------------------------------------------
function [zeros, err, fZ] = findzeros(xA, step, xE, fname)
   numZeros = 0;             % Anzahl bisher gefundener Nullstellen.
   A = NaN;                  % Matrix mit Intervallen um Nullstellen.
   X = (xA:step:xE)';        % Intervallgrenzen festlegen.
   if X(rows(X) < xE)        % Intervall endet mit xE!
      X(rows(X)+1) = xE;
   endif
   Y = arrayfun(fname, X);   % Funktionswerte bestimmen.
   i = 1;                    % Vorzeichenwechsel finden.
   LL = rows(X);
   SGN = sign(Y(1));
   while (i <= LL)
      if (SGN != sign(Y(i))) % Vorzeichenwechsel! Nur für i > 1 möglich.
         SGN = sign(Y(i));
         A(++numZeros,1) = X(i-1);
         A(  numZeros,2) = X(i);
      endif
      i++;
   endwhile
   if (numZeros < 1)         % Keine Vorzeichenwechsel gefunden.
      error "No zeros found. Try smaller x_step."
   endif
   % Matrix A enthält nun je Zeile ein Intervall, in dem sich zumindest
   % eine Nullstelle befindet. Wir finden je Intervall eine Nullstelle
   % per Bisektion.
   i = 1;
   LL = rows(A);
   while (i <= LL)
      a = A(i,1);  b = A(i,2);
      F(i) = feval(fname, A(i,1));
      F_b = feval(fname, A(i,2));
      while true
         c = (a + b) / 2;
         F_c = feval(fname, c);
         if ((c == a) || (c == b))
            break;
         elseif (F_c == 0)
            a = c;  b = c;
            break;
         elseif (sign(F_c) == sign(F(i)))
            a = c;
         else
            b = c;
         endif
      endwhile
      A(i,1) = a;  A(i,2) = b;
      F(i) = F_c;
      i++;
   endwhile
   zeros = A(:,1);
   err = (A(:,2) - A(:,1));
   fZ  = F;
endfunction