function [ t, x, u, l, cost ] = search13 ( a, b, A, B, d1, d2, M, x0, T, N )

%*****************************************************************************80
%
%% search13() performs optimal control for case #13.
%
%  Discussion:
%
%    Predator-Prey Model
%
%    minimize
%      N1(T) + integral ( 0 <= t <= T ) A / 2 u(t)^2 dt
%    subject to
%      N1'(t) = ( 1 - N2(t) ) N1(t) - d1 N1(t) u(t)
%      N2'(t) = ( N1(t) - 1 ) N2(t) - d2 N2(t) u(t)
%      z'(t) = u(t)
%      N1(0) = N10
%      N2(0) = N20
%      z(0) = 0
%      z(T) = B
%   with Hamiltonian:
%      H = A / 2 u(t)^2 
%        + lambda1(t) ( ( 1 - N2(t) ) N1(t) - d1 N1(t) u(t) ) 
%        + lambda2(t) ( ( N1(t) - 1 ) N2(t) - d2 N2(t) u(t) )
%        + lambda3(t) u(t)
%    with adjoint equation:
%      lambda1(t)' = - lambda1(t) ( 1 - N2(t) - d1 u(t) ) - lambda2  N2(t)
%      lambda1(T) = 0.0
%      lambda2(t)' = + lambda1(t) N1(t) + lambda2(t) ( N1(t) - 1 - d2 u(t) )
%      lambda2(T) = 0.0
%      lambda3(t)' = ?
%      lambda3(T) = input
%    and optimality condition:
%      A u(t) - d1 N1(t) lambda1(t) - d2 N2(t) lambda2(t) + lambda3(t)
%    and
%      0 <= u(t) <= M
%      integral ( 0 <= t <= T ) u(t) dt = B
%
%  Modified:
%
%    10 October 2022
%
%  Author:
%
%    Suzanne Lenhart, John Workman
%    Modifications by John Burkardt
%
%  Reference:
%
%    Suzanne Lenhart, John Workman,
%    Optimal Control Applied to Biological Models,
%    Chapman and Hall, 2007
%    LC: QH324.8.L46,
%    ISBN: 978-1-58488-640-2
%
%  Input:
%
%    real a, b: two starting values for theta.
%
%    real A: a weight coefficient for the cost function.
%
%    real B: the total pesticide application.
%
%    real d1, d2: death rate coefficients.
%
%    real M: the pesticide application level.
%
%    real x0(3): the initial values for the state variables.
%
%    real T: the final time.
%
%    integer N: the number of time steps.
%
%  Output:
%
%    real t(N+1): the time values.
%
%    real N1(N+1), N2(N+1): the final state estimate.
%
%    real u(N+1), z(N+1): the final optimal control estimate.
%
%    real l(N+1): the final adjoint estimate.
%
%    real cost(*): the cost history.
%
  flag = -1;

  [ t, x, u, l, cost ] = opt13 ( a, A, d1, d2, M, x0, T, N );
  Ta = x(N+1,3) - B;
  disp('1 iteration completed.');

  [ t, x, u, l, cost ] = opt13 ( b, A, d1, d2, M, x0, T, N );
  Tb = x(N+1,3) - B;
  disp('2 iterations completed.')

  n = 2;

  while ( n < 100 )

    n = n + 1;

    if ( abs ( Tb ) < abs ( Ta ) )
      t = a;
      a = b;
      b = t;
      t = Ta;
      Ta = Tb;
      Tb = t;
    end

    d = Ta * ( b - a ) / ( Tb - Ta );
    b = a;
    Tb = Ta;
    a = a - d;
    [ t, x, u, l, cost ] = opt13 ( a, A, d1, d2, M, x0, T, N );
    Ta = x(N+1,3) - B;
    
    if ( abs ( Ta ) < 1E-10 )
      flag = 1;
      break;
    end

    n = n + 1;
    fprintf ( 1, '  %d  %g  %g\n', n, a, Ta );

  end

  fprintf ( 1, '  Final estimate was a = %g\n', a );

  return
end

