//
// Prima parte
//
function [z,s,d,iter] = Bisezione(f,a,b)
    //
    // Applica il metodo di bisezione alla funzione f a partire dall'intervallo
    // [a,b]. L'iterazione si arresta quando l'intervallo non varia oppure
    // f(z) = 0.
    //
    // z: approssimazione suggerita.
    // s,d: estremo sinistro e destro, rispettivamente, dell'ultimo
    //      intervallo calcolato.
    // iter: numero di iterazioni effettuate.
    //
    iter = 0;
    if f(a)*f(b) > 0 then
       error('la funzione non assume valori di segno opposto agli estremi');
    else
        z = (a+b)/2;
        while ( z <> a & z <> b & f(z) <> 0 ),
              if sign(f(a)) <> sign(f(z)) then b = z;
              else a = z; end;
              z = (a+b)/2;
              iter = iter + 1;
        end;
        s = a; d = b;
    end;
endfunction
//
function y = F(x)
    y = (x - 2).^13;
endfunction
//
function y = G(x)
    p = poly(2*ones(1,13),'x');
    y = horner(p,x);
endfunction
//
clf(); 
a = 0;
b = 3;
[zf,s,d,iter] = Bisezione(F,a,b)
//
printf('\n   zf = %10.9e , F(zf) = %3.2e , iter = %d , ampiezza = %3.2e\n', ...
       zf,F(zf),iter,d-s);
halt('...premere INVIO per continuare (1)...');
a = 0; //+ rand();
b = 3; //+ rand();
[zg,s,d,iter] = Bisezione(G,a,b)
//
printf('\n   zg = %10.9e , G(zg) = %3.2e , iter = %d , ampiezza = %3.2e\n', ...
       zg,G(zg),iter,d-s);
halt('...premere INVIO per continuare (2)...');
//
x = linspace(-1,1,901)'; x = x + 2;
plot2d(x,[G(x),F(x)],style = [5,2]); xgrid();
xlabel('x'); legend('G(x)','F(x)','in_upper_left');
halt('...premere INVIO per continuare (3)...');
//
clf();
x = linspace(-0.23,0.23,700)'; x = x + 2;
plot2d(x,[G(x),F(x)],style = [5,2]); xgrid();
xlabel('x'); legend('G(x)','F(x)','in_upper_left');
assi = gca();
assi.children(2).children(1).thickness=3;
plot2d(zg,G(zg),style = -5);
pause;
//
// Seconda parte
//
function x = MetodoNewton(f,d1f,x0,delta)
    kmax = 500;
    k = 0;
    x = x0;
    while (abs( f(x($)) / d1f(x($)) ) >= delta & k < kmax),
        x($+1) = x($) - f(x($)) / d1f(x($));
        k = k + 1;
    end;
endfunction
//
function y = d1F(x)
    y = 13 * (x - 2).^12;
endfunction
//
clf();
x0 = 3;
delta = 1d-5;
zN = MetodoNewton(F,d1F,x0,delta);
printf('\n   zN = %10.9e , F(zN) = %3.2e , iter = %1.0f , zN - 2 = %3.2e\n', ...
       zN($),F(zN($)),length(zN)-1, zN($) - 2);
plot2d(x,F(x),style = 5);
xlabel('x'); ylabel('F(x)'); xgrid();
plot(zN($),F( zN($) ), 'ko');
//
//  Esercizio:
//
//  Adattare quanto fatto nella prima parte al caso:
//
//   f(x) = (3*x - 2)^13
//
