function r = newton(f,fp,x0,abstol,maxits)
 x = x0;
 for i = 1:maxits
    dx = - f(x)/fp(x);
    x = x + dx;
    fprintf('x(%i)=%20.16e\n',i,x)
    if abs(dx) < abstol
        r = x;
        fprintf('Converged to desired accuracy in %i steps',i)
        return;
    end
 end
 error('Failed to converge in %i steps',maxits)