function [x f conv] = newtfun(fh, dfh, x0)
% newtfun solves f(x)=0 using Newton's Methond
% Input arguments
%   fh    handle to f(x)
%   dfh   handle to f'(x)
%   x0    initial guess
% Output arguments
%   x     root
%   f     f(x)
%   conv  convergence=1, 0 otherwise

format long e
steps = 0; x = x0; re = 1e-8;

myrel = 1;

while (myrel > re) & (steps < 20)
    xold = x;
    x = x - fh(x)/dfh(x);
    steps = steps + 1;
    disp([x fh(x)])
    myrel = abs((x-xold)/x);
end

if myrel <= re
    conv = 1;
else
    conv = 0;
end

f = fh(x);