function x = fixpoint(g, guess, eboundrel, dispflag)
%fixpoint - find a fixed point of the given function g.
%  if "dispflag" is nonzero, print a summary on each iteration
%   The value returned, x, is either:
%      a real value, such that (g(x) - x)/g(x) < eboundrel, or
%      NaN, if we diverge
%   We use two checks for divergence:
%      - values getting further apart on successive iterations
%      - the effective magnitude of the slope >1 on successive iterations
DIVMAX = 3;
x = guess;
prev = Inf;
prevdiff = Inf;
slopecount = 0;  % number of iterations with slope > 1
diffcount = 0;   % number of iterations getting further apart
itercount = 0;   % number of iterations, for testing
ea = Inf;
if dispflag
    fprintf('Iter        x           g(x)           ea.         diff      slope\n');
    fprintf('----   ------------  ------------   ---------      ----      -----\n');
end
while ea > eboundrel && slopecount < DIVMAX && diffcount < DIVMAX
    prev2 = prev; % save two back for slope computation
    prev = x;    % save the current value
    x = g(x);   % compute the next one
    diff = abs(x-prev);
    slope = abs(diff/(prev-prev2));  %x==g(prev) & prev==g(prev2)
    % check for divergence
    if slope >= 1
       slopecount = slopecount + 1;
    else
       slopecount = 0;
    end
    % check for divergence
    if diff >= prevdiff  % if steps are not getting closer together 
       diffcount = diffcount + 1;  % too many times in a row
    else
        diffcount = 0;      % reset whenever it gets smaller
    end
    if dispflag
        fprintf('%4d  %12.8f  %12.8f  %12.8f  %10.6f %6.2f\n',itercount,prev,x,ea,diff,slope);
    end
    prevdiff = diff;
    ea = abs(diff/x);
    itercount = itercount + 1;
end
if ea > eboundrel  % alas, we diverged...
    x = NaN;
end
end

