function bvp_shooting_test ( )

%*****************************************************************************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.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    18 August 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Josef Stoer, Roland Bulirsch,
%    Introduction to Numerical Mathematics,
%    Springer, 2002,
%    ISBN: 038795452X,
%    LC: QA297.S8213.
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'bvp_shooting_test\n' );
  fprintf ( 1, '  A simple example of the shooting method,\n' )
  fprintf ( 1, '  in which a boundary value problem (BVP) is solved\n' );
  fprintf ( 1, '  using an initial value problem (IVP) solver,\n' );
  fprintf ( 1, '  a guess for the missing initial condition ALPHA\n' );
  fprintf ( 1, '  and a nonlinear equation solver.\n' );
%
%  The BVP is specified over [a,b]:
%
  a = 0.0;
  b = 1.0;
%
%  The side conditions are y(a) = ya, y(b) = yb.
%
  ya = 4.0;
  yb = 1.0;
%
%  Define F(alpha) as the computed value of y(b) - yb.
%  Evaluate F at two starting points, then take 10 secant steps.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Step  Alpha  YB(Alpha)  F(ALPHA)\n' );
  fprintf ( 1, '\n' );

  for i = -1 : 10

    if ( i == -1 )
      alpha = 0.0;
    elseif ( i == 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

    [ x, y ] = ode45 ( @myode, [a;b], [ya; alpha ] );
    yb_alpha = y(end,1);
    f_alpha = yb_alpha - yb;
    fprintf ( '  %d  %g  %g  %g\n', i, alpha, yb_alpha, f_alpha );

  end
%
%  Now plot the solution.
%
  sol = ode45 ( @myode, [a,b], [ya, alpha ] );
  x = linspace ( a, b, 101 );
  y = deval ( sol, x );
  plot ( x, y(1,:), 'Linewidth', 3 );
  grid ( 'on' );
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y(X) --->' );
  title ( 'Solution of y"=3/2 y^2, y(0)=4, y(1)=1.' )
  filename = 'bvp_shooting_test.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );
  
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'bvp_shooting_test\n' );
  fprintf ( 1, '  Normal end of execution.\n' )
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end
function dydx = myode ( x, y )

%*****************************************************************************80
%
%% myode evaluates the right hand side of the differential equation.
%
%  Discussion:
%
%    y" = 3/2 y^2 is rewritten as a pair of first order equations:
%
%    y'(1) = y(2)
%    y'(2) = 3/2 y(1)^2
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    18 August 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real x(*): values of the independent variable.
%
%    real y(2,*): values of the dependent variable.
%
%    real dydx(2,*): values of the derivative.
%
  dydx(1,:) = y(2,:);
  dydx(2,:) = 1.5 * y(1,:).^2;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% timestamp prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end

