function [ t, x, u, lambda, cost ] = test05 ( A, B, C, x0 )

%*****************************************************************************80
%
%% test05()
%
%  Discussion:
%
%    Maximize
%      integral ( 0 <= t <= 1 ) A x - B u^2 dt
%    subject to
%      x' = -1/2 x^2 + C u
%    x(0) = x0 > -1, A, B > 0.
%
%  Modified:
%
%    30 September 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real A, B, C: three problem parameters.
%
%    real X0: the prescribed value of X at time 0.
%
%  Output:
%
%    real T(*), X(*), U(*), LAMBDA(*), COST: the time, state, control,
%    adjoint, and cost.
%
  delta = 0.001;
  n = 1001;
  t0 = 0.0;
  t1 = 1.0;

  step = 0;
  step_max = 50;

  while ( step <= step_max )

    if ( step == 0 )

      t = linspace ( t0, t1, n );
      x = zeros(1,n);
      u = zeros(1,n);
      lambda = zeros(1,n);
      cost = - Inf;

    else
%
%  Forward Euler applied to dxdt = - 1/2 x^2 + C u
%
      x_old = x;
      dt = ( t1 - t0 ) / ( n - 1 );
      for i = 1 : n
        if ( i == 1 )
          x(i) = x0;
        else
          x(i) = x(i-1) + dt * ( -0.5 * x(i-1)^2 + C * u(i-1) );
        end
      end
%
%  Forward Euler applied (in backwards time) to d lambda dt = -A + lambda x
%
      lambda_old = lambda;
      lambda0 = 0.0;
      dt = ( t0 - t1 ) / ( n - 1 );
      for i = n : -1 : 1
        if ( i == n )
          lambda(i) = lambda0;
        else
          lambda(i) = lambda(i+1) + dt * ( - A + lambda(i+1) * x(i+1) );
        end
      end
%
%  Update the control
%  Use lambda to estimate U, but average with previous result.
%
      u_old = u;
      u1 = C * lambda / ( 2.0 * B );
      u = 0.5 * ( u1 + u_old );
%
%  Evaluate COST.
%
      cost_old = cost;
      f = A * x + B * u.^2;
      cost = trapz ( t, f );
%
%  Convergence tests.
%
      disp ( [ step, norm(u), norm(x), norm(lambda), abs(cost) ] );

      temp1 = norm ( u - u_old ) / norm ( u );
      temp2 = norm ( x - x_old ) / norm ( x );
      temp3 = norm ( lambda - lambda_old ) / norm ( lambda );
      temp4 = abs ( j - j_old ) / abs ( j );
      test = min ( temp1, min ( temp2, min ( temp3, temp4 ) ) );
      disp ( [ step, temp1, temp2, temp3, temp4 ] );
      if ( 10 < step && test < delta )
        return;
      end

    end

    step = step + 1;

  end

  fprintf ( 1, '\n' );
  fprintf ( 1, 'test05(): Warning!\n' );
  fprintf ( 1, '  Maximum number of steps taken without convergence.\n' );

  return
end

