function [ a, b, x, job ] = bisection_rc ( a, b, fx, job )

%*****************************************************************************80
%
%% bisection_rc() seeks a zero of f(x) in a change of sign interval.
%
%  Discussion:
%
%    The bisection method is used.
%
%    This routine uses reverse communication, so that the function is always
%    evaluated in the calling program.
%
%    On the first call, the user sets JOB = 0, and the values of A and B.
%    Thereafter, the user checks the returned value of JOB and follows 
%    directions.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 January 2013
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A, B, the endpoints of the change of sign interval.  
%    The user only initializes these values on first call.
%
%    real FX, the function value at X.
%    Set this to 0 before the first call.  On subsequent calls, FX
%    should be set to the function value at X.
%
%    integer JOB, a communication flag.
%    Set JOB to 0 before the first call.
%
%  Output:
%
%    real A, B, the updated endpoints of the change of sign interval.  
%
%    real X, a point at which the function is to 
%    be evaluated.
%
%    integer JOB, a communication flag.
%    The program sets this to 1.
%    
  persistent fa
  persistent fb
  persistent state
  persistent x_local
%
%  JOB==0: Initialization.
%
  if ( job == 0 )

    fa = 0.0;
    fb = 0.0;
    state = 1;
    x = a;
    job = 1;
%
%  STATE==1: 
%    Received function value at X=A.
%    Request function value F(B).
%
  elseif ( state == 1 )

    fa = fx;
    x = b;
    state = 2;
%
%  STATE == 2.
%    Received function value at X=B.
%    Request function value at midpoint.
%  
  elseif ( state == 2 )

    fb = fx;

    if ( sign ( fa ) == sign ( fb ) )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'BISECTION_RC - Fatal error!\n' );
      fprintf ( 1, '  F(A) and F(B) have the same sign.\n' );
      error ( 'BISECTION_RC - Fatal error!' );
    end

    x = ( a + b ) / 2.0;
    state = 3;
%
%  STATE == 3.
%    Received function value at X=B.
%    Request function value at midpoint.
%  
  else

    if ( sign ( fx ) == sign ( fa ) )
      a = x_local;
      fa = fx;
    else
      b = x_local;
      fb = fx;
    end
    x = ( a + b ) / 2.0;
    state = 3;

  end
%
%  Remember X.
%
  x_local = x;

  return
end

