function [x,nit,xv,flag]=newton(fun,dfun,x)
nmax=1000;
toll=1.e-12;
flag=1;
xv=[x];
for nit=1:nmax
    inc=feval(fun,x)/feval(dfun,x);
    x=x-inc;
    xv=[xv x];
    if abs(inc)<toll, flag=0; return, end
end
