// METODO DELLE SUCCESSIVE BISEZIONI 
// PER LA SOLUZIONE DI EQUAZIONI NON LINEARI
//
// function [x,fx,exitflag,output] = bzero(fun,interv,options)
//
// DATI DI INPUT:
//   funz = stringa contenente il nome della funzione
//   interv = [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-1   numero di iterate maggiore di nmax
//                         oppure inf o nan
//        
//    output = struttura
//       output.algorithm = algoritmo usato
//        output.funcCount  = valutazioni di funzioni
//        output.iterations = numero di iterazioni
//
function [x,fx,exitflag,output,xvet] = bzero(fun,interv,options)

if length(interv) ~= 2
    error(' interv deve avere due elementi' )
end 
   
output.algorithm = 'bisezioni';
a = interv(1);
b = interv(2);
fa = feval(a,fun);
fb = feval(b,fun);
if abs(fa) == 0
   x = a;
   fx = fa;
   exitflag =1;
   output.funcCount   = 2;
   output.iterations= 0;
    return
end
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
nargin=argn(2);
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 = (options.Display);
            select 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


fc = 1/%eps;
nit = 1;
fcont=2;
if displayit
    mprintf('Val. funzioni Iterationi            x           f(x)\n')
end
xvet = [];

while ( abs(b-a)/max(1,min(abs(a),abs(b)))>tol ) & (nit <= nmax)
    c = (a+b)*0.5;
    xvet = [xvet c];
    fc = feval(c,fun);
    fcont = fcont+1;
      if displayit
       mprintf('%10.0f %10.0f   %13.6g %13.6g    \n',fcont, nit,c,fc')
    end
    if abs(fc) ~= 0
        
    if (fc*fa < 0)
        b = c;
        fb = fc;
    else
        a = c;
        fa = fc;
    end
  
    nit = nit+1;
   else
       
        x = c;
       fx = fc;
   exitflag =1;
   output.funcCount   = fcont;
   output.iterations= nit;
  
   return
   end
end

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