function [ t, x, u, l, cost ] = opt14 ( B, k, x0, N )

%*****************************************************************************80
%
%% opt14() performs optimal control for case #14.
%
%  Discussion:
%
%    Invasive Plant Species
%
%    minimize
%      sum ( 1 <= j <= 5 ) r(j,N+1) + B sum ( 1 <= t <= N ) u(j,t)^2 )
%    subject to the difference equation:
%      r(j,t+1) = ( r(j,t) + k r(j,t) / ( epsilon + r(j,t) ) ( 1 - u(j,t) )
%      r(j,1) = rho(j)
%    with Hamiltonian:
%      H = ?
%    with adjoint equation:
%      ?
%    and optimality condition:
%      ?
%    and
%      0 <= u(j,t) <= 1
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  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 B: a weight coefficient in the cost function.
%
%    real k: the spread rate.
%
%    real x0(5): initial values for the state variables.
%
%    real N: the number of time steps.
%
%  Output:
%
%    real t(N+1): the times, from 0 to T, in unit steps.
%
%    real x(N+1,5): the state variables.
%
%    real u(N,5): the control variables.
%
%    real l(N+1,5): the adjoint variables.
%
%    real cost(): the cost history.
%
  t0 = 0.0;
  t1 = N;
  epsilon = 0.01;
%
%  Initialize the problem variables.
%
  t = linspace ( t0, t1, N+1 );
  t = t';
  tb = t(end:-1:1);
  x = zeros(N+1,5);
  x(1,1:5) = x0;
  u = zeros(N+1,5);
  l = zeros (N+1,5);
  l1 = [ 1.0, 1.0, 1.0, 1.0, 1.0 ];
  l(N+1,1:5) = l1;
  cost = [];
%
%  Forward/backward sweep.
%
  it = 0;

  while ( true )
%
%  Save previous variables for convergence test.
%
    oldu = u;
    oldx = x;
    oldl = l;
%
%  Update state variables.
%
    for i = 1 : N
      x(i+1,1:5) = ( x(i,1:5) + k * x(i,1:5) ./ ( epsilon + x(i,1:5) ) ) ...
        .* ( 1.0 - u(i,1:5) );
    end
%
%  Update adjoint variables.
%
    for i = N : -1 : 1
      l(i,1:5) = l(i+1,1:5) .* ( 1.0 - u(i,1:5) ) ...
        .* ( 1.0 + k * epsilon ./ ( epsilon + x(i,1:5) ).^2 );
    end
%
%  Update the controls.
%
    u(1:N,1:5) = 0.5 / B * l(2:N+1,1:5) ...
      .* ( x(1:N,1:5) + k * x(1:N,1:5) ./ ( epsilon + x(1:N,1:5) ) );   

    u = max ( u, 0.0 );
    u = min ( u, 1.0 );
    u = 0.5 * ( u + oldu );
%
%  Evaluate the cost.
%
    fint = sum ( x(N+1,1:5) ) + B * sum ( sum ( u(1:N,1:5).^2 ) );
    cost = [ cost, fint ];
%
%  Convergence test.
%
    ratio = Inf;
    ratio = convergence ( ratio, x(:,1), oldx(:,1) );
    ratio = convergence ( ratio, x(:,2), oldx(:,2) );
    ratio = convergence ( ratio, x(:,3), oldx(:,3) );
    ratio = convergence ( ratio, x(:,4), oldx(:,4) );
    ratio = convergence ( ratio, x(:,5), oldx(:,5) );
    ratio = convergence ( ratio, u(:,1), oldu(:,1) );
    ratio = convergence ( ratio, u(:,2), oldu(:,2) );
    ratio = convergence ( ratio, u(:,3), oldu(:,3) );
    ratio = convergence ( ratio, u(:,4), oldu(:,4) );
    ratio = convergence ( ratio, u(:,5), oldu(:,5) );
    ratio = convergence ( ratio, l(:,1), oldl(:,1) );
    ratio = convergence ( ratio, l(:,2), oldl(:,2) );
    ratio = convergence ( ratio, l(:,3), oldl(:,3) );
    ratio = convergence ( ratio, l(:,4), oldl(:,4) );
    ratio = convergence ( ratio, l(:,5), oldl(:,5) );
    it = it + 1;
    if ( 0.0 <= ratio )
      fprintf ( 1, 'opt14(): Convergence after %d iterations\n', it );
      break;
    end

  end

  return
end

