function [ x, y ] = bvp_shooting ( a, b, ya, yb, myode )

%*****************************************************************************80
%
%% bvp_shooting_test() shows the shooting method for boundary value problems.
%
%  Discussion:
%
%    Demonstrate how the solution of a boundary value problem (BVP) can be
%    approximated using an initial value problem (IVP) solver, using the
%    shooting method.
%
%    The idea is to guess the missing initial condition ALPHA, use the IVP
%    solver to get the corresponding solution, call F(ALPHA) the difference
%    between the computed solution and the desired end condition, and then
%    employ a nonlinear equation solver to seek ALPHA such that F(ALPHA)=0.
%
%    The secant iteration is repeated up to ITER_MAX times, a value which
%    can be adjusted by the user.
%
%    The secant iteration is terminated early if the desired boundary
%    condition is approximated to a tolerance of TOL, a value which can
%    be adjusted by the user.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 March 2026
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Josef Stoer, Roland Bulirsch,
%    Introduction to Numerical Mathematics,
%    Springer, 2002,
%    ISBN: 038795452X,
%    LC: QA297.S8213.
%
%  Input:
%
%    real a, b: the left and right endpoints of the interval.
%
%    real ya, yb: the left and right boundary conditions:
%    y(a) = ya, y(b) = yb.
%
%    function myode: evaluates the right hand side of the BVP.
%
%  Output:
%
%    real x(*), y(*): the vectors of data defining the computed solution.
%
  iter_max = 10;
  tol = 1.0E-06;

  tspan = [ a, b ];
%
%  Define F(alpha) as the computed value of y(b) - yb.
%  Evaluate F at two starting points, then take secant steps.
%  We hope that we can compute a value of alpha such that the
%  initial value problem with y'(a)=alpha will solve the boundary
%  value problem with y(b) = yb.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Iteration  Alpha  YB(Alpha)  F(ALPHA)\n' );
  fprintf ( 1, '\n' );

  for iter = -1 : iter_max

    if ( iter == -1 )
      alpha = 0.0;
    elseif ( iter == 0 )
      beta = alpha;
      f_beta = f_alpha;
      alpha = -1.0;
    else
      gamma = alpha - f_alpha * ( beta - alpha ) / ( f_beta - f_alpha );
      beta = alpha;
      f_beta = f_alpha;
      alpha = gamma;
    end

    y0 = [ ya; alpha ];
%
%  In this call to ode45(), we don't use an AT sign on the function, 
%  because it's already been "AT-ed" in the calling function...
%  I doubt that statement is clear the first time you read it.
%
    [ x, y ] = ode45 ( myode, tspan, y0 );
    yb_alpha = y(end,1);
    f_alpha = yb_alpha - yb;
    fprintf ( '  %d  %g  %g  %g\n', iter, alpha, yb_alpha, f_alpha );

    if ( abs ( f_alpha ) < tol )
      break;
    end

  end

  return
end
