function [ a, b, fa, fb ] = nonlin_bisect ( f_string, a, b )

%*****************************************************************************80
%
%% nonlin_bisect() is an interactive bisection function.
%
%  Discussion:
%
%    This function manages an interactive computation which uses bisection
%    to seek a zero of a function f(x) in the interval a <= x <= b.
%
%    The user enters a formula for f(x), and the values of a and b.
%
%    The program can be invoked by a function call, in which case the
%    string specifying f(x) must be quoted:
%      [ a, b, fa, fb ] = nonlin_bisect ( 'x^2-4', -1, 3 )
%    or, if called with no arguments, it will request them:
%      [ a, b, fa, fb ] = nonlin_bisect ( );
%      Enter function formula, like x^2-4: x^2-4
%      Enter left limit, a: -1
%      Enter right limit, b: 3
%
%    The function is specified as a string which is either:
%    * a MATLAB expression using the argument 'x', 
%    * the name of an M-file followed by the argument '(x)'.
%
%    The string should not contain any spaces between symbols, except
%    when it is passed as a function argument in quotes.
%
%    It is not necessary to use the "dot" notation for expressions involving
%    '*', '/', or '^', but it doesn't hurt either.
%
%    Examples of function specifications:
%
%      x^2
%      x.^2
%      3/(x^4+5*x-6)
%      sin(7*x)*sqrt(x)/8
%      wiggle(x)     <-- where "wiggle.m" is a user-provided M file.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    25 July 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_string: the function whose zero is desired.
%
%    real a, b: the change of sign interval.
%
%  Output:
%
%    real a, b: the reduced interval.
%
%    real fa, fb: the values f(a) and f(b).
%

%
%  Get F, A, B.
%
  if ( nargin < 1 )
    f_string = input ( 'Enter function formula, like x^2-4: ', 's' );
  end
  f_element = vector_to_element ( f_string );
  f = str2func ( [ '@(x)', f_element ] );

  if ( nargin < 2 )
    a = input ( 'Enter left limit a: ' );
  elseif ( ischar ( a ) )
    a = str2num ( a );
  end

  if ( nargin < 3 )
    b = input ( 'Enter right limit b: ' );
  elseif ( ischar ( b ) )
    b = str2num ( b );
  end
%
%  Refuse if no change of sign.
%
  if ( sign ( f(a) ) == sign ( f(b) ) )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'nonlin_bisect - Fatal error!\n' );
    fprintf ( 1, '  f(a) and f(b) have the same sign.\n' );
    error ( 'nonlin_bisect - Fatal error!' );
  end
%
%  Carry out bisection.
%
  tol = 10.0 * eps * ( b - a );
  [ a, b, it ] = bisect ( a, b, tol, f );
  x = ( a + b ) / 2.0;
  fa = f(a);
  fb = f(b);
  fx = f(x);

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Function = "%s"\n', f_string );
  fprintf ( 1, '  a = %24.16g, f(a) = %g\n', a, fa );
  fprintf ( 1, '  b = %24.16g, f(b) = %g\n', b, fb );
  fprintf ( 1, '  Interval tolerance = %g\n', tol );
  fprintf ( 1, '  Number of bisections = %d\n', it );
  fprintf ( 1, '  x = %24.16g, f(x) = %g\n', x, fx );

  return
end
function [ a, b, it ] = bisect ( a, b, tol, f )

%*****************************************************************************80
%
%% bisect() carries out the bisection method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    17 June 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A, B, the endpoints of a change of sign interval.
%
%    real TOL, the interval size tolerance.  Once |B-A| < TOL,
%    the iteration will stop.
%
%    pointer value=F(x), a pointer to a function.  This can be
%    an anonymous function or the name of a MATLAB M-file function.
%
%  Output:
%
%    real A, B, the new endpoints that constitute an
%    change of sign interval no larger than TOL.
%
%    integer IT, the number of bisections.
%
  fa = f ( a );
  fb = f ( b );
  it = 0;

  while ( tol < abs ( b - a ) )
    c = ( a + b ) / 2.0;
    fc = f(c);
    it = it + 1;
    if ( sign ( fc ) == sign ( fa ) )
      a = c;
      fa = fc;
    else
      b = c;
      fb = fc;
    end
  end

  return
end
function s2 = s_escape_tex ( s1 )

%*****************************************************************************80
%
%% s_escape_tex() de-escapes TeX escape sequences.
%
%  Discussion:
%
%    In particular, every occurrence of the characters '\', '_',
%    '^', '{' and '}' will be replaced by '\\', '\_', '\^',
%    '\{' and '\}'.  A TeX interpreter, on seeing these character
%    strings, is then likely to return the original characters.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    09 February 2010
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string S1, the string to be de-escaped.
%
%  Output:
%
%    string S2, a copy of the string, modified to avoid TeX escapes.
%
  s1_length = length ( s1 );

  s1_pos = 0;
  s2_pos = 0;
  s2 = [];

  while ( s1_pos < s1_length )

    s1_pos = s1_pos + 1;

    if ( s1(s1_pos) == '\' || ...
         s1(s1_pos) == '_' || ...
         s1(s1_pos) == '^' || ...
         s1(s1_pos) == '{' || ...
         s1(s1_pos) == '}' )
      s2_pos = s2_pos + 1;
      s2 = strcat ( s2, '\' );
    end

    s2_pos = s2_pos + 1;
    s2 = strcat ( s2, s1(s1_pos) );

  end

  return
end
function f_element = vector_to_element ( f_vector )

%*****************************************************************************80
%
%% vector_to_element() rewrites a MATLAB vector function elementwise.
%
%  Discussion:
%
%    In MATLAB, by default, the operators '*', '/', and '^' are interpreted
%    as linear algebraic operations.  In many cases, a user simply desired
%    the elementwise operation instead, which MATLAB indicates using a 
%    dotted notation: '.*', './', and '.^'.
%
%    This function converts a string to a form that uses dotted operators.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 June 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_vector, a string presumably representing a MATLAB expression.
%
%  Output:
%
%    string f_element, a copy of f_vector in which all operators have been
%    replaced by their dotted (elementwise) form.
%
  n_vector = length ( f_vector );
  f_element = [''];

  c_prev = '';
  for iv = 1 : n_vector
    c = f_vector(iv);
    if ( ( c == '*' || c == '/' || c == '^' ) && ( c_prev ~= '.' ) )
      f_element = [ f_element, ' .', c, ' ' ];
    else
      f_element = [ f_element, c ];
    end
    c_prev = c;
  end

  return
end

