function [Ms,M,Md]=maxabsfun(a,b,fun,dfun,trgr,Ng,varargin)

% détermination du max de la valeur absolue d'une fonction réelle
%
% *********************************************************
%  [Ms,M,Md]=maxabsfun(a,b,fun,dfun,trgr,Ng,varargin) :
%
%  calcul de M=max |f(x)| pour x décrivant [a,b]
%
%  * Variables d'entrée : 
%     *  a,b : bornes de l'intervalle d'étude (ils peuvent être
%              symboliques)
%     * fun  : fonction f (vectorizée)
%  * Variables d'entrée optionnelles (peuvent être vides ou inexistantes) :
%     * dfun   : fonction f' (si absente est calculée en symbolique)
%     * trgr   : entier dans {0,1} avec graphique si =1 (=0 par défaut)
%     * Ng     : entier pour le nombre de points pour le graphiques et le
%               calcul du maximum discret. égal à 1000 par défaut
%     * varargin : arguments optionnels de la fonction f (et f')
%  * Variables de sortie : 
%     * Ms : maximum M en symbolique (est vide si la recherche symbolique
%            n'a pas fonctionné)
%     * M  : maximum M
%     * Md : maximum discret, c'est-à-dire le max de f sur
%            x=linspace(a,b,Ng)
%
% ************ Fonctions auxiliaires utilisées ************
%
%    annee_version 
%
% *********************************************************
%
% EXEMPLES
%
% [Ms,M,Md]=maxabsfun(0,2*pi,'sin',[],1);
% [Ms,M,Md]=maxabsfun(0,2*pi,inline('sin(a*x)','x','a'),inline('a*cos(a*x)','x','a'),1,[],1);
% [Ms,M,Md]=maxabsfun(-.2,1,inline('x.^4+x.^3-48*x.^2+56*x+12'),[],1);
% [Ms,M,Md]=maxabsfun(0,1,inline('exp(-x.^2).*(16*x.^4-48*x.^2+12)'),[],1);
% [Ms,M,Md]=maxabsfun(0,2,inline('exp(-x.^2).*(16*x.^4-48*x.^2+12)'),[],1);
% [Ms,M,Md]=maxabsfun(2,4,inline('-6^(1/4).*(x-5)./(((-x.^2+10*x-6).^(3/4)))/2'),[],1);
%
% syms x;
% f=diff(exp(-x.^2),4);
% f=inline(vectorize(char(f)));
% [Ms,M,Md]=maxabsfun(0,1,f,[],1);
%
%
% 2014 by          Jérôme BASTIEN
%                  Université Claude Bernard Lyon I, UFRSTAPS, Laboratoire CRIS, Villeurbanne
%                  E-Mail : jerome.bastien@univ-lyon1.fr

nin=nargin;
nout=nargout;
if nin<=3
    testdf=1;
else
    testdf=isempty(dfun);
end
if nin<=4
    trgr=0;
else
    if isempty(trgr)
        trgr=0;
    end
end
if nin<=5
    Ng=1000;
else
    if isempty(Ng)
        Ng=1000;
    end
end
syms x;
if testdf
    dfuns=diff(feval(fun,x,varargin{:}),x);
else
    dfuns=feval(dfun,x,varargin{:});
end
%S=unique(solve(char(dfuns)));
if annee_version>=2019
    S=unique(solve(dfuns));
else
    S=unique(solve(char(dfuns)));
end
if ~isempty(S)
    Sd=double(S);
    indd=abs(imag(Sd))<=eps&real(Sd)>double(a)&real(Sd)<double(b);
    S=S(indd);
    if ~isempty(S)
        if trgr
            xg=double(real(S));
            Mg=feval(fun,xg,varargin{:});
        end
        F=[S.',sym(a),sym(b)];
        MM=abs(feval(fun,F,varargin{:}));
        [M,ind]=max(double(MM));
        ind=ind(1);
        Ms=MM(ind);
    end
end
if isempty(S)
    options=optimset('TolX',eps);
    funbis=@(x) feval(fun,x,varargin{:});
    xm=fminbnd(funbis,double(a),double(b),options);
    funbis=@(x) -feval(fun,x,varargin{:});
    xM=fminbnd(funbis,double(a),double(b),options);
    F=[xm,xM,a,b];
    if trgr
        xg=[xm,xM];
        Mg=feval(fun,xg,varargin{:});
    end
    MM=double(abs(feval(fun,F,varargin{:})));
    [M,ind]=max(MM);
    ind=ind(1);
    M=MM(ind);
    Ms=[];
end
if trgr
    x=linspace(double(a),double(b),Ng);
    y=feval(fun,x,varargin{:});
    clf;
    mmm=min(y);
    MMM=max([y,abs(y)]);
    ec=MMM-mmm;
    hold on;
    plot(x,y,x,abs(y),xg,Mg,'r*',[double(a),double(b)],[M,M],'k');
    legend('fonction','valeur absolue','extrémaux','max','location','best');
    ylim([mmm,MMM+.1*ec]);
    if mmm<0
        plot([double(a),double(b)],[0,0],'-.k');
    end
    hold off;
    if nout==3
        Md=max(abs(y));
    end
else
    if nout==3
        x=linspace(double(a),double(b),Ng);
        y=feval(fun,x,varargin{:});
        Md=max(abs(y));
    end
end