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