function fx = p01_fun ( option, nvar, x )

%*****************************************************************************80
%
%% p01_fun() evaluates the function for problem 1.
%
%  Title:
%
%    The Freudenstein-Roth function
%
%  Description:
%
%    One way to use a continuation code as a nonlinear root finder
%    is to start with a set of nonlinear equations G(X), and an
%    approximate root A, and create a "homotopy" function F(X,Y)
%    with the properties that F(A,0.0) = 0 and F(X,1.0) = G(X).
%    Thus, the homotopy function F has a known exact solution
%    from which we can start with no difficulty.  If the continuation
%    code can take us from Y = 0 to Y = 1, then we have found
%    an X so that F(X,1.0) = 0, so we have found a solution to G(X)=0.
%
%    The Freudenstein-Roth function F(X) is derived in this way
%    from a homotopy of G(X):
%
%      F ( X(1), X(2), X(3) ) =
%        G ( X(1), X(2) ) - ( 1 - X(3) ) * G ( Y1, Y2 )
%
%    where Y1 and Y2 are some fixed values, and
%
%      G(1) = X(1) - X(2)*X(2)*X(2) + 5*X(2)*X(2) -  2*X(2) - 13
%      G(2) = X(1) + X(2)*X(2)*X(2) +   X(2)*X(2) - 14*X(2) - 29
%
%  Options 1, 2, 3:
%
%    The starting point is X0 = ( 15, -2, 0 ).
%
%    A great deal of information is available about the homotopy curve
%    generated by this starting point:
%
%    The function F(X) has the form
%
%      F(1) = X(1) - X(2)**3 + 5*X(2)**2 -  2*X(2) - 13 + 34*(X(3)-1)
%      F(2) = X(1) + X(2)**3 +   X(2)**2 - 14*X(2) - 29 + 10*(X(3)-1)
%
%    There is a closed form representation of the curve in terms of the
%    second parameter:
%
%      X(1) = (-11*X(2)**3 + 4*X(2)**2 + 114*X(2) + 214) /  6
%      X(2) = X(2)
%      X(3) = (    X(2)**3 - 2*X(2)**2 -   6*X(2) +   4) / 12
%
%    The first option simply requests the production of solution points
%    along the curve until a point is reached whose third component is
%    exactly 1.
%
%    Options 2 and 3 use the same starting point, and also stop when the
%    third component is 1.  However, these options in addition search
%    for limit points in the first and third components of the solution,
%    respectively.
%
%    The target solution has X(3) = 1, and is ( 5, 4, 1 ).
%
%    Limit points for X1:
%
%      ( 14.28309, -1.741377,  0.2585779 )
%      ( 61.66936,  1.983801, -0.6638797 )
%
%    Limit points for X3:
%
%     (20.48586, -0.8968053, 0.5875873)
%     (61.02031,  2.230139, -0.6863528)
%
%    The curve has several dramatic bends.
%
%
%  Options 4, 5, and 6:
%
%    The starting point is (4, 3, 0).
%
%    The function F(X) has the form
%
%      F(1) = X(1) - X(2)**3 + 5*X(2)**2 -  2*X(2) - 13 +  3*(X(3)-1)
%      F(2) = X(1) + X(2)**3 +   X(2)**2 - 14*X(2) - 29 - 31*(X(3)-1)
%
%    There is a closed form representation of the curve in terms of the
%    second parameter:
%
%      X(1) = (14*X(2)**3 -79*X(2)**2 +52*X(2) + 245) / 17
%      X(2) = X(2)
%      X(3) = (   X(2)**3 - 2*X(2)**2 - 6*X(2) +   9) / 17
%
%    The correct value of the solution at X(3)=1 is:
%
%      (5, 4, 1)
%
%    In option 5, limit points in the first component are sought,
%    and in option 6, limit points in the third component are
%    sought.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    02 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Ferdinand Freudenstein, Bernhard Roth,
%    Numerical Solutions of Nonlinear Equations,
%    Journal of the Association for Computing Machinery,
%    Volume 10, 1963, Pages 550-556.
%
%  Input:
%
%    integer OPTION, the option index.
%
%    integer NVAR, the number of variables.
%
%    real X(NVAR), the argument of the function.
%
%  Output:
%
%    real FX(NVAR-1), the value of the function at X.
%

%  Get the starting point, Y.
%
  y = p01_start ( option, nvar );
%
%  G is the function value at the starting point,
%  F the function value at the current point.
%
  gy = p01_gx ( y );

  gx = p01_gx ( x );
%
%  The parameter X3 generates the homotopy curve.
%
  fx(1:nvar-1) = gx(1:nvar-1) + ( x(3) - 1.0 ) * gy(1:nvar-1);

  return
end
