function [xvect,xdif,fx,nit]=newton(x0,toll,nmax,fun,dfun,trgr,varargin)

% [xvect,xdif,fx,nit]=newton(x0,toll,nmax,fun,dfun,trgr,varargin)
%
% Newton method for nonlinear systems
%
%  >> INPUT <<
%
%  x0         starting point
%  toll       tolerance required
%  nmax       maximum number of iterations
%  fun        function given (must be vectorized !)
%
%  >> OPTIONNAL INPUT <<
%
%  dfun       first derivate of this function
%             dfun can be empty or non-existent
%             In this case, dfun is determinated thanks to the symbolic
%             calculus
%  trgr       integer : 0 : no graph is drawn
%                     : 1 : a graph with the different values is drawn
%             trgr can be empty or non-existent; in this case, no graph is
%             drawn
%  varargin   common additional variables of functions fun and dfun
%             Warning, roots of function fun with respect to first
%             variable of function fun
%
%
%  >> OUTPUT <<
%
%  xvect      vector containing all the iterates, the root is
%             xvect(end)
%  xdif       difference between to iterates
%  fx         values of fun in xvect
%  nit        number of iterations required to reach the
%             prescribed tolerance
%
%  >> EXAMPLE
%  [xvect,xdif,fx,nit]=newton(0.7,1e-8,100,inline('sin(2*x)-1+x'),inline('2*cos(2*x) + 1'),1)
%  [xvect,xdif,fx,nit]=newton(0.7,1e-8,100,inline('sin(2*x)-1+x'))
%  [xvect,xdif,fx,nit]=newton(0.7,1e-8,100,inline('sin(2*x)-1+x'),[],1)
%  [xvect,xdif,fx,nit]=newton(0.7,1e-8,5,inline('x.^2-a','x','a'),[],[],2)
%  [xvect,xdif,fx,nit]=newton(0.7,1e-8,5,inline('x.^2-a','x','a'),[],1,2)
%  [xvect,xdif,fx,nit]=newton(0.7,1e-8000,5000,inline('x.^2'),[],1)

% ligne extraite de quad, par exemple (est-ce pertinent dans la mesure où
% elle ne vectorize pas toujours !)
fun=fcnchk(fun);
nin=nargin;
if nin<=4
    testdf=1;
else
    testdf=isempty(dfun);
end
if nin<=5
    trgr=0;
else
    if isempty(trgr)
        trgr=0;
    end
end
if testdf
    syms x;
    g=diff(feval(fun,x,varargin{:}),x);
    nb=length(varargin);
    % AP : 
    % Améliorer cela TROP LENT dans le cas vardf=0
    if nb==0
        dfun=inline(char(g));
        dfun=fcnchk(dfun);
        vardf=1;
    else
        vardf=0;
    end
else
    vardf=1;
end
err=toll+1;
nit=0;
xvect=zeros(1,nmax+1);
fx=xvect;
xdif=zeros(1,nmax);
xvect(1)=x0;
fx(1)=feval(fun,x0,varargin{:});
x=x0;
while (nit < nmax & err > toll & fx(nit+1)~=0)
    nit=nit+1;
    % AP : 
    % Améliorer cela TROP LENT dans le cas vardf=0
    if vardf
        dfx=feval(dfun,x,varargin{:});
    else
        dfx=g;
        if ~(isnumeric(dfx));
            dfx=eval(dfx);
        end
    end
    if (dfx == 0)
        err=toll*1.e-10;
        disp('Stop for vanishing dfun');
    else
        xn=x-fx(nit)/dfx;
        err=abs(xn-x);
        xdif(nit)=err;
        x=xn;
        xvect(nit+1)=x;
        fx(nit+1)=feval(fun,x,varargin{:});
    end;
end;
xvect=xvect(1:nit+1);
fx=fx(1:nit+1);
xdif=xdif(1:nit);
if trgr
    ng=1000;
    xxmin=min(xvect);
    xxmax=max(xvect);
    if xxmin==xxmax
        xxmin=x0-1;
        xxmax=x0+1;
    end
    ecasr=.1*(xxmax-xxmin);
    xxmin=xxmin-ecasr;
    xxmax=xxmax+ecasr;
    if nit~=0
        X=xvect(1:nit);
        Xp=X;
        X=X(ones(1,2),:);
        X=(X(:)).';
        X=[X,x(end)];
        Yp=zeros(1,nit);
        Y=[Yp;fx(1:nit)];
        Y=(Y(:)).';
        Y=[Y,0];
    else
        X=x(end);
        Y=0;
    end
    xg=sort([linspace(xxmin,xxmax,ng),xvect]);
    %         AP : 
    %         Essayer de vectorizer en utilisant, qui doit valable pour 'sin' et @sin
    %         nb=length(varargin);
    %         if nb==0
    %         fung=inline(vectorize(char(fun)));
    %         yg=feval(fung,xg,varargin{:});
    %         ou essayer aussi avec exist et whos
    %         FAIRE la même chose dans la fonction fixepoint
    yg=feval(fun,xg,varargin{:});
    clf;
    hold on;
    if nit~=0
        plot(xg,yg,'b',Xp,Yp,'og',X(end),Y(end),'*r');
        legend('fonction','itérés de Newton','dernière valeur',0);
    else
        plot(xg,yg,'b',X(end),Y(end),'*r');
        legend('fonction','unique valeur',0);
    end
    plot(X,Y,'r',[xxmin,xxmax],[0,0],'k:');
    xlim([xxmin,xxmax]);
    hold off;
end