function [x, fx, a, b] = newtroot(fname, a, b, xtol, ftol, fig, start)
% Newton's method for a differentiable function F mapping the real line to 
% the real line, with bisection safeguards
% input
%   fname : string giving name of function that evaluates F(x) and
%      its derivative F'(x), with call: [f, deriv] = fname(x)
%   a, b: points satisfying a<b and F(a), F(b) have opposite signs
%   xtol: tolerance on length of interval containing root
%   ftol: tolerance on f (default 1e-8)
%   fig: figure number for plot, or 0 for no plot (default 0)
%   start: -1 for starting at left end, 1 for starting at right end
% output
%   x: best point found
%   fx: F(x)
%   a: left end point of final interval containing root of F
%   b: right end point of final interval containing root of F
if nargin < 4
    xtol = 1e-8;
end
if nargin < 5
    ftol = 1e-8;
end
if nargin < 6
    fig = 0;
end
if nargin < 7
    start = -1;
end
if a >= b
    error('newtroot: a >= b')
end
[fa, deriva] = feval(fname, a);
[fb, derivb] = feval(fname, b);
if fa*fb > 0
    error('newtroot: F(a), F(b) have same sign')
end
if start < 0
    x = a; 
    fx = fa;
    deriv = deriva;
else
    x = b;
    fx = fb;
    deriv = derivb;
end
if fig > 0
    figure(fig)
    hold on
    plot(a, fa, 'r+')
    plot(b, fb, 'r+')
end
maxit = 20;
iter = 0;
while b-a > xtol & abs(fx) > ftol & iter < maxit
    iter = iter + 1;
    newtonstep = -fx/deriv;  % OK to divide by zero
    % make sure newton step doesn't go the wrong way or go past
    % the bisection step
    if (x == a & newtonstep < 0) | (x == b & newtonstep > 0) |...
     abs(newtonstep) > (b-a)/2
        x = (a + b)/2;
        fprintf('bisection step to %g', x)
    else
        x = x + newtonstep;
        fprintf('Newton step to %20.15e', x);
    end
    [fx, deriv] = feval(fname, x);
    if fx*fa > 0 % update a 
        a = x;
        fa = fx;
    else % update b
        b = x;
        fb = fx;
    end
    fprintf(', new values a= %g fa = %g b = %g fb = %g\n', a, fa, b, fb)
    if fig > 0
        plot(x, fx, 'm*')
    end
end
    
    
    
