// Graphe

deff('[y]=myF(x)','y=x.^5-2*x.^3-8*x.^2+16');

clf();
subplot(2,2,1)
x=-1.75:0.1:2.5;
plot(x,myF(x),x,zeros(x));
title('Graphe de F');

// Détermination des zéros

format(20) // grande précision

x1=fsolve(-1.5,myF,10^{-20})
Fprime1=mygradF(x1)

x2=fsolve(1,myF,10^{-20})
Fprime2=mygradF(x2)

x3=fsolve(3,myF,10^{-20})
Fprime3=mygradF(x3)

pause;

// Newton

deff('[y]=mygradF(x)','y=5*x.^4-6*x.^2-16*x');

function X=myNewton(x0,N)
    X=[x0];
    y=x0;
    for i=1:N
        y=y-mygradF(y)\myF(y);
        X=[X y];
    end
endfunction;

// Tracé des approximations

subplot(2,2,2)

plot(0:20,x1*ones(0:20));
plot(0:20,x2*ones(0:20));
plot(0:20,x3*ones(0:20));

title('Suite des approximations');

plot(0:20,myNewton(-5,20),'ro',0:20,myNewton(1,20),'k*',0:20,myNewton(5,20),'g.')

pause;

// Ordre de convergence
// NB: le nombre d'itérations a été arrêtée avant que l'erreur ne soit trop petite et déclenche une erreur de la fonction log 

X=myNewton(-5,9);
Y=log(abs(X-x1));
Z=Y(2:10)./Y(1:9);

subplot(2,2,3)
plot(2:10,Z,'ro');
m=mean(Z(6:9))
plot2d(2:10,m*ones((2:10)),rect=[2,-5,10,4],style=[1]);
title('Rapport du log des erreurs partant de -5');

subplot(2,2,4)
[a b]=reglin(6:10,log(Y(6:10)));
plot(1:10,log(Y),'ro',1:10,a*(1:10)+b,'b');
title('Log du log des erreurs partant de -5');
exp(a)

printf('Les ordres de convergence numériques sont '+string(m)+' et '+string(exp(a)))