function [aa,bb,wr,flag]=bisezione(f,a,b,toll,maxit)

% Algoritmo di bisezione, con criterio di arresto sul residuo pesato e
% ampiezza dell'intervallo.
%
% Dati di ingresso:
% f: funzione (inline function)
% a: estremo sinistro
% b: estremo destro
% toll: tolleranza richiesta per il test del residuo pesato
% maxit: massimo indice dell'iterata permesso
%
% Dati di uscita:
% aa: sequenza degli estremi di sinistra degli intervalli [a_k,b_k], 
%     immagazzinata in vettore colonna;
% bb: sequenza degli estremi di destra degli intervalli [a_k,b_k],
%     immagazzinata in vettore colonna;
% wr: sequenza dei residui pesati,immagazzinata in vettore colonna;
% flag: 0 processo terminato correttamente, 
%     1 processo non terminato correttamente.

if b < a, s=b; b=a; a=s; end % Aggiusta errori utente.

flag=0;  fa=feval(f,a); fb=feval(f,b); 
aa=[a]; bb=[b]; wr=[];

% test: ipotesi bisezione non soddisfatte
if fa*fb > 0, flag=1; return; end

% zero all'estremo iniziale "a".
if fa == 0, aa=[a]; bb=[a]; return; end % a sol.

% zero all'estremo iniziale "b".
if fb == 0, aa=[b]; bb=[b]; return; end % b sol.



% iterazioni bisezione
for k=1:maxit
    c=(a+b)/2;  fc=feval(f,c); % punto medio di [a_k,b_k]
    w=(b-a)/(fb-fa); % peso "w".
    wres=abs(fc*w);  % residuo pesato.
    wr=[wr; wres];   % aggiorna sequenza residui pesati
    % uscita per raggiungimento risultato.
    if (wres<toll) |(fc==0) 
        aa(k+1)=c; bb(k+1)=c; return;% OK exit.
    end 
    % determinazione intervallo [a_{k+1},b_{k+1}]
    if fc*fa > 0 % "c" sostituisce "a"
        % aggiorna "aa", "bb", "a", "fa"
        aa(k+1)=c; bb(k+1)=b; a=c; fa=fc;
    else % "c" sostituisce "b"
        % aggiorna "aa", "bb", "b", "fb"
        aa(k+1)=a; bb(k+1)=c; b=c; fb=fc;
    end 
end

% se si e' raggiunto questo punto, si sono fatte troppe iterazioni
flag=1;
