function [r,fr]=newton(fh,dfh,x0)
%inputs:fh=f(x),dfh=f'(x), x0=initial guess
%outputs:x=approx root,f=f(x)
steps=0;
x=x0;
re=1e-8;
myrel=1;
while myrel>re && (steps<20)
    steps=steps+1;
    xold=x;
    x=x-fh(x)/dfh(x);
    myrel=abs((x-xold)/x);
end
r=x;
fr=fh(r);
end