function ZEROS = muller ( F, FNREAL, ZEROS, NPREV, MAXIT, EP1, EP2)

%*****************************************************************************80
%
%% muller() implements Muller's method.
%
%  pages 122-4 of conte/deboor
%
%  I N P U T
%
%  F  the handle of the function whose  N  zeros are being sought
%  FNREAL  a logical variable. If true, all approximations are taken
%     to be real. allowing this routine to be used even if  F(Z)  is
%     defined only for real  Z .
%  ZEROS(1:NPREV)  previously found zeros (if any)
%  ZEROS(NPREV+1:END) first guesses for the remaining zeros to be
%     found. (If you know nothing, then 0 is as good a guess as any.)
%  MAXIT  maximum number of function evaluations allowed per zero.
%  EP1  iteration is stopped if abs(H) < EP1*abs(ZR), with
%     H = latest change in zero estimate  zero .
%  EP2  although the  EP1  criterion is not met, iteration is stopped if
%     abs(F(zero)) < EP2.
%  NPREV   number of zeros previously found.
%
% O U T P U T
%
%  ZEROS   the vector of zeros found
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
               % initialization
   EPS1 = max(EP1, 1e-12); EPS2 = max(EP2, 1e-20);

for I=NPREV+1:length(ZEROS)
   KOUNT = 0;
        % compute the first three estimates for zero as
        %     ZEROS(i)+.5, ZEROS(i)-.5, ZEROS(i)
   IFLAG = true;
   while IFLAG 
      ZERO = ZEROS(I);
      [FZR,DVDF1P,KOUNT,ZEROS,IFLAG] = dflate (F, ZERO+.5, I, KOUNT, ZEROS);
      % The fortran program uses a control structure (indicated by the
      % input argument *[number]) which permits specification inside the
      % subroutine that, on leaving, the next statement to be executed be
      % the one in the calling program labeled [number] (rather than the
      % one following the calling statement).
      % Since MATLAB eschews go to statements, this made it impossible to
      % copy more or less the Fortran coding here.
      if ~IFLAG, IFLAG = true; continue, end
      [FZR,FZRPRV,KOUNT,ZEROS,IFLAG] = dflate (F, ZERO-.5, I, KOUNT, ZEROS);
      if ~IFLAG, IFLAG = true; continue, end
      HPREV = -1; DVDF1P = (FZRPRV - DVDF1P)/HPREV;
      [FZR,FZRDFL,KOUNT,ZEROS,IFLAG] = dflate(F, ZERO, I, KOUNT, ZEROS); 
      if ~IFLAG, IFLAG = true; continue, end
      H = .5; 
      while KOUNT<=MAXIT  
         %       or  H  is not small
         %       or  FZR = F(ZERO)  is not small
         %       or  FZRDFL = FdflateD(ZERO)  is not small or not much
         %           bigger than its previous value  FZRPRV
         %  the following is done
         DIVDF1 = (FZRDFL - FZRPRV)/H; DIVDF2 = (DIVDF1 - DVDF1P)/(H+HPREV);
         HPREV = H; DVDF1P = DIVDF1; 
         C = DIVDF1 + H*DIVDF2;
         SQR = C*C - 4*FZRDFL*DIVDF2; if FNREAL && real(SQR)<0, SQR  = 0; end
         SQR = sqrt(SQR);
         if real(C)*real(SQR) + imag(C)*imag(SQR) < 0, DEN = C - SQR;
         else                                          DEN = C + SQR; end
         if abs(DEN) <= 0, DEN = 1; end
         H = -2*FZRDFL/DEN; FZRPRV = FZRDFL; ZERO = ZERO + H;
         if KOUNT > MAXIT, IFLAG = false; continue, end
         goon = true;
         while goon
            [FZR,FZRDFL,KOUNT,ZEROS,IFLAG] = dflate(F, ZERO, I, KOUNT, ZEROS);
            if ~IFLAG, IFLAG = true; goon = false; continue, end
                % check for convergence
            if abs(H) < EPS1*max(abs(ZERO),1),
                IFLAG = false; goon = false;  KOUNT=MAXIT+1; continue, end
            if max(abs(FZR), abs(FZRDFL)) <= EPS2,
                IFLAG = false; goon = false;  KOUNT=MAXIT+1; continue, end
                % check for divergence
            if abs(FZRDFL) < 10*abs(FZRPRV), goon = false; continue, end
            H = H/2; ZERO = ZERO - H;
         end
      end
      ZEROS(I) = ZERO;
   end
end

function [FZERO,FZRDFL,KOUNT,ZEROS,IFLAG] = dflate( F, ZERO, I, KOUNT, ZEROS)

KOUNT = KOUNT+1; FZERO = F(ZERO); FZRDFL = FZERO;IFLAG = true;
if I<2, return, end
for j=2:I
   DEN = ZERO - ZEROS(j-1);
   if abs(DEN) == 0, ZEROS(I) = ZERO*1.001; IFLAG = false; return, end
   FZRDFL = FZRDFL/DEN;
end
return
