function [rac,nb_iter]=muller_elem(init,exp_f,arret)
% Produit une racine d'équation (E) f(x)=0 par méthode de Muller.
% Le calcul se termine par un forçage à une écriture symbolique de rac,
% qui est affichée en plus de la numérique!
% Ceci peut etre aisément modifié!
%
% variables d'entrée
% init est un vecteur de trois valeurs initiales distinctes
% pour la méthode de Muller  classique. On peut généraliser!
% exp_f est l'expression de la fonction f sous la forme 'f(x)';
% arret est un vecteur de trois réels [f_seuil,nul_seuil,n_max]
% dont la signification est la suivante:
% si |f(x)|<f_seuil on considère f(x) nul;
% si |x|<nul_seuil on considère x nul;
% n_max est un entier naturel qui désigne le nombre maximal
% d'itérations autorisées.
%
% variables de sortie
% rac est une valeur approchée de la racine obtenue à partir d'init
% par algorithme de Muller classique; elle peut etre réelle ou complexe.
% nb_iter désigne le nombre d'itérations exécutée pour sortir rac.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Fonctions connexes appelées
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spécification
format long e;

% Tests d'entrée

if nargin~=3
    error('Nombre des arguments incorrect, ou écriture inattendue');
end
% Les tests de controle suivants peuvent etre plus ou moins développés,
% selon la dangerosité pressentie de l'utilisateur... ordinaire!
if (size(init,1)~=1)|(size(init,2)~=3)|((sum(diff(init)==0)>0))
    error('Champ init mal écrit. Voir controle');
end
if (size(arret,1)~=1)|(size(arret,2)~=3)
    error('Champ arret mal écrit. Voir controle');
end
% compléments éventuels


% Algorithme proprement dit

% Etape 1:initialisations
X=init;Y=zeros(1,3);

%  La suite serait avantageusement remplacée par une écriture vectorielle
% mais l'utilisateur calculera faux sans le voir..; Il faudrait écrire
% 'x.^2+1' au lieu de 'x^2+1'
for k=1:3
    x=X(k);Y(k)=eval(exp_f);
end

d=diff(X);dif_div=diff(Y)./d;
nb_iter=0;
bool_arret=0;

% Etape 2
while bool_arret==0
    nb_iter=nb_iter+1;
    % détermination du polynome d'interpolation de f sur le support x
    % et résolution de l'équation approchée
    An=diff(dif_div)/sum(d);
    Bn=dif_div(2)+An*d(2);
    x=X(3);Cn=eval(exp_f);
    if abs(Cn)<arret(1)
        bool_arret=1;
    else
        % détermination du bon epsilon fournissant le dénominateur optimal
        den=Bn+sqrt(Bn^2-4*An*Cn).*[1 -1];[val,ind]=max(abs(den));
        den_opt=den(ind);
        ajout=-2*Cn/den_opt;
        
        % préparation de la boucle suivante
        
        X(1:2)=X(2:3);X(3)=X(2)+ajout;
        Y(1:2)=Y(2:3);
        x=X(3);Y(3)=eval(exp_f);
        d=diff(X);dif_div=diff(Y)./d;
        %mise à jour de bool_arret
        bool_arret=...
            (nb_iter>arret(3))|(abs(ajout)<arret(2))|(abs(Y(3))<arret(1));
       
    end
end

% toilettage de la racine:
% on la force réelle si elle l'est presque (au sens de nul_seuil)
% si elle est complexe, on renvoie aussi sa conjuguée dont on sait
% qu'elle est aussi solution par théorie des corps si l'expression
% de f est à coefficient réels.

if abs(real(X(3)))<arret(2)
    X(3)=imag(X(3))*i;
end
if abs(imag(X(3)))<arret(2)
    rac=real(X(3));
else
    rac=[X(3) conj(X(3))];
end

% sortie des résultats sous forme symbolique
% Libérer du commentaire, si souhaité, les deux lignes suivantes.
% Voir aussi la fonction affiche_racines.

disp('Sortie des résultats sous forme symbolique proposée par matlab');
symb=sym(rac,'r');disp('rac = ');disp (symb);

% fin de fonction