//
function y = f(x)
    y = x.*(x - 1).*(x - 2) - c
endfunction
//
function dy = d1f(x)
    dy = 3*x.^2 - 6*x + 2
endfunction
//
exec('/home/ciampa/Lavoro/Scilab/Lavoro/Zeri/Metodi/newton1d.sci');
E_newt = 1d-6;
kmax = 100;
//
clf();
N = 100;
cmax = 1;
c = 0;
[z, v, info] = newton1d(f,1d-3,d1f,E_newt,kmax,'zitto');
if info ~= 1 then printf('\n  c = %3.2e   ;   info = %d',c,info);
                  plot(z,c,'ro');
                  else plot(z,c,'go'); end;
for k = 1:N,
    c = k*cmax/N;
    [z, v, info] = newton1d(f,z,d1f,E_newt,kmax,'zitto');
    if info ~= 1 then printf('\n  c = %3.2e   ;   info = %d',c,info);
                      plot(z,c,'ro');
                  else plot(z,c,'go'); end; 
end;
c = cmax;
w = z;
for k = 1:N,
    c = (N-k)*cmax/N;
    [w, v, info] = newton1d(f,w,d1f,E_newt,kmax,'zitto');
    if info ~= 1 then printf('\n  c = %3.2e   ;   info = %d',c,info);
                      plot(w,c,'r+');
                  else plot(w,c,'b+'); end;
end;
xgrid();
xlabel('z');
ylabel('c');
titolo = 'kmax = ' + string(kmax);
xtitle(titolo);
