function test03 ( )

%*****************************************************************************80
%
%% test03() seeks an optimal control of a model problem.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 September 2022
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'test03():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Seek an optimal control for a model biological problem.()\n' );

  t0 = 0.0;
  t1 = 5.0;
  n = 100;
  t = linspace ( t0, t1, n + 1 );
  u = zeros ( n + 1, 1 );
  x = zeros ( n + 1, 2 );
  lam = zeros ( n + 1, 2 );

  mu = + 0.1;
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Control function update rate mu = %g\n', mu );

  for ufun = 2 : 2

    fprintf ( 1, '\n' );
    fprintf ( 1, '  Using control function U #%d\n', ufun );

    for step = 1 : 40

%     fprintf ( 1, '    Step #%d\n', step );
%
%  On first step, control function 1 or 2 is specified.
%  On subsequent steps, current control function is updated.
%
      if ( step == 1 )
        if ( ufun == 1 )
          u(1:n+1) = 0.0 * ( 0.0 <= t(1:n+1) & t(1:n+1) <  1.0 ) ...
            + 0.5 * ( 1.0 <= t(1:n+1) & t(1:n+1) <  3.0 ) ...
            + 1.0 * ( 3.0 <= t(1:n+1) & t(1:n+1) <= 5.0 );
        elseif ( ufun == 2 )
          u(1:n+1) = 1.0 * ( 0.0 <= t(1:n+1) & t(1:n+1) <  2.0 ) ...
            + 0.5 * ( 2.0 <= t(1:n+1) & t(1:n+1) <  4.0 ) ...
            + 0.0 * ( 4.0 <= t(1:n+1) & t(1:n+1) <= 5.0 );
        end
      else
        u(1:n+1) = u(1:n+1) ...
          + mu * ( lam(1:n+1,1) - lam(1:n+1,2) ) .* x(1:n+1,1);
        u = max ( u, 0.0 );
        u = min ( u, 1.0 );
      end
%
%  Forward Euler on x.
%
      dt = ( t1 - t0 ) / n;
      for i = 1 : n + 1
        if ( i == 1 )
          x(i,1:2) = [ 1.0, 2.0 ];
        else
          x(i,1) = x(i-1,1) + dt *         u(i-1)   * x(i-1,1);
          x(i,2) = x(i-1,2) + dt * ( 1.0 - u(i-1) ) * x(i-1,1);
        end
      end
%
%  Evaluate cost.
%
      cost = trapz ( t, log(x(:,2)) );
      fprintf ( '    Step %d, Cost = %g\n', step, cost );
%
%  Integrate backward in time for LAMBDA.
%
      dt = ( t0 - t1 ) / n;
      for i = n + 1 : -1 : 1
        if ( i == n + 1 )
          lam(i,1:2) = [ 0.0, 0.0 ];
        else
          lam(i,1) = lam(i+1,1) ...
            + dt * ( - lam(i+1,1) * u(i+1) - lam(i+1,2) * ( 1.0 - u(i+1) ) );
          lam(i,2) = lam(i+1,2) ...
            + dt * ( - 1.0 / x(i+1,2) );
        end
      end
%
%  Display the current control.
%
      plot ( t, u, 'linewidth', 3 );
      grid ( 'on' );
      xlabel ( '<-- t -->' );
      ylabel ( '<-- u(t) -->' );
      s = sprintf ( 'Control #%d, step %d', ufun, step );
      title ( s );
      filename = sprintf ( 'test03/u%d_step%d.png', ufun, step );
      print ( '-dpng', filename );

    end

    plot ( t, x, 'linewidth', 3 );
    grid ( 'on' );
    xlabel ( '<-- t -->' );
    ylabel ( '<-- X(t) -->' );
    s = sprintf ( 'State variable X(t) for control #%d step #%d', ufun, step );
    title ( s );
    legend ( '{x_1}', '{x_2}' );
    filename = sprintf ( 'test03/x_u%d_step%d.png', ufun, step );
    print ( '-dpng', filename );

    plot ( t, lam, 'linewidth', 3 );
    grid ( 'on' );
    xlabel ( '<-- t -->' );
    ylabel ( '<-- {\lambda}(t) -->' );
    s = sprintf ( 'Adjoint {\\lambda}(t) for control #%d, step %d', ufun, step );
    title ( s );
    legend ( '{\lambda_1}', '{\lambda_2}' );
    filename = sprintf ( 'test03/lam_u%d_step%d.png', ufun, step );
    print ( '-dpng', filename );

  end

  return
end

