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

%--------------------------------------------------------------------------
% Oggetto: 
% metodo della "regula falsi" per risolvere equazioni nonlineari "f(x)=0".
%--------------------------------------------------------------------------
% Input:
% f: funzione continua di cui si vuole calcolare uno zero.
% a, b: estremi dell'intervallo iniziale [a,b] comprendente lo zero, 
%       tale che "f(a1)*f(b1) < 0"
% toll: tolleranza richiesta
% maxit: numero massimo di iterazioni.
%--------------------------------------------------------------------------
% Output:
% aa,bb: vettori la cui k-sima componente contiene "a(k)" e "b(k)", 
%      tali che "f(a(k))*f(b(k)) < 0",
% cc: valori "c(k)" proposti dal metodo di regula-falsi.
% wr: vettore dei residui pesati.
% flag: 0 terminazione corretta, 1 terminazione per troppe iterazioni.
%--------------------------------------------------------------------------

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=[]; cc=[];

% test: ipotesi regula falsi 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 regula falsi
for k=1:maxit
    
    % Nota: a questa altezza sono noti "a", "b", "fa", "fb".
    c=(a*fb-b*fa)/(fb-fa); % iterata regula falsi
    cc=[cc; c]; % update vettore dei vari "c"
    
    fc=feval(f,c); % valutazione "f(c)"
    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;


