% METODO DELLE SUCCESSIVE BISEZIONI E SECANTI
% PER LA SOLUZIONE DI EQUAZIONI NON LINEARI
%
% function [x,fx,exitflag,output] = bzero(fun,int,options)
%
% DATI DI INPUT:
%   funz = stringa contenente il nome della funzione
%   int = [a,b]  estremi dell'intervallo in cui si cerca la radice di funz
%   options = struttura contenente dati opzionali 
%     options.TolX     = tolleranza
%     options.Nmax    = numero massimo di iterate
%     options.Display = livello di visualizzazione 
%                'off'    nessuna visualizzazione;
%                'iter'   visualizzazione per ogni iterata; 
%                'final'  visualizzazione finale;
%                'notify' (default) visualizzazione solo se il metodo non converge.
% DATI DI OUTPUT:
%    x   = soluzione
%    fx  = valore della funzione in x
%    exitflag = flag di errore 
%         exitflag = 1  nessun errore
%         exitflag = -1  numero di iterate maggiore di nmax
%        
%    output = struttura
%       output.algorithm = algoritmo usato
%        output.funcCount  = valutazioni di funzioni
%        output.iterations = numero di iterazioni
%
function [x,fx,exitflag,output] = bzero(fun,int,options)

if length(int) ~= 2
    error(' Int deve avere due elementi' )
end 
   
output.algorithm = 'bisezioni, secanti';
a = int(1);
b = int(2);
fa = feval(fun,a);
if abs(fa) == 0
   x = a;
   fx = fa;
   exitflag =1;
   output.funcCount   = 1;
   output.iterations= 0;
    return
end
fb = feval(fun,b);
if abs(fb) == 0
   x = b;
   fx = fb;
   exitflag =1;
   output.funcCount   = 2;
   output.iterations= 0;
    return
end

if fa*fb > 0 
    error('La funzione non cambia di segno agli estremi')
end
if nargin == 2 
    tol  = eps; 
    nmax = 100;
    displaynot = 1;
    displayit  =0;
    displayfin =0;
else
    if isfield(options,'TolX')
        if ~isempty(options.TolX)
            tol=options.TolX;
        else
            tol=eps;
        end
    else
        tol = eps;
    end
    if isfield(options,'Nmax')
        if ~isempty(options.Nmax)
            nmax = options.Nmax;
        else
            nmax=100;
        end
    else
        nmax=100;
    end
    if isfield(options,'Display')
        if ~isempty(options.Display)
            display = lower(options.Display);
            switch display
            case {'off'}
                displaynot =0;
                displayit  =0;
                displayfin =0;   
            case{'final'}
                displaynot =0;
                displayit  =0;
                displayfin =1;
            case{'iter'}
                displaynot =0;
                displayit  =1;
                displayfin =0;
            case{'notify'}
                displaynot =1;
                displayit  =0;
                displayfin =0;
            end
        else
            displaynot =1;
            displayit  =0;
            displayfin =0;
        end
    else
        displaynot =1;
        displayit  =0;
        displayfin =0;
    end
end


nit = 1;
fcont=2;
if displayit
    fprintf('Val. funzioni Iterationi            x           f(x)\n')
end
x0 = a;
fx0 = fa;
x1 = b;
fx1 = fb;

while ( abs(x0-x1)/max(1,min(abs(x0),abs(x1)))>tol ) & (nit <= nmax)
  
    c = x1 - fx1*(x1-x0)/(fx1-fx0);
    fc = feval(fun,c);
    fcont = fcont+1;
  
 
      
   
    
      if  c < a |  c > b
        c = a + (b-a)/2;
        fc = feval(fun,c);
        fcont = fcont+1;
        if fa*fc < 0
         b = c; fb = fc;
         x1 = c; fx1 = fc;
         x0 = a; fx0 = fa;
        else
         a=c; fa=fc;
         x1 = c; fx1=fc;
         x0 = b; fx0=fb;
        end
        metodo = 'B';
    else
         if fa*fc < 0
         b = c; fb = fc;
         x0 = x1; fx0=fx1;
         x1 = c; fx1 = fc;
        else
         a=c; fa=fc;
         x0 = x1; fx0=fx1;
         x1 = c; fx1=fc;
        end
        metodo = 'S';
    end
     if abs(fc) == 0 
       x = c;
       fx = fc;
       exitflag =1;
       output.funcCount   = fcont;
       output.iterations= nit;
       return
    end
      if displayit
          fprintf('%10.0f %10.0f   %13.6g %13.6g  %s  \n',fcont, nit,c,fc,metodo')
      end
    nit = nit+1;
   
  
end

if nit > nmax | isnan(c) | isinf(c) | isnan(fc) | isinf(fc)
   x=c;
   fx=fc;
    exitflag =-1;
   output.funcCount   = fcont;
   output.iterations= nit;    
   if displaynot 
       if nit > nmax
        fprintf(' Numero massimo di iterate raggiunto \n')
    else
         fprintf(' Inf o NaN')
    end
    
        fprintf('Val. funzioni Iterationi            x           f(x)\n') 
        fprintf('%10.0f %10.0f   %13.6g %13.6g    \n',fcont, nit,x,fx')   
    end
else
   x = c;
   fx = fc;
   exitflag =1;
   output.funcCount   = fcont;
   output.iterations= nit;
   if displayfin
      fprintf('Soluzione trovata nell''intervallo [%13.6g %13.6g]',a,b)
       fprintf('Val. funzioni Iterationi            x           f(x)\n')
      fprintf('%10.0f %10.0f   %13.6g %13.6g    \n',fcont, nit,x,fx')
    end
   
end
