function y = shoot( f, ya, yb, z1, z2, t, tol )
  % --------------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % Written: April 15, 1999
  % Slightly Revised: April 17, 2019
  %
  % $Id: shoot.m,v 1.5 2007/04/18 15:14:57 senning Exp $
  %
  % Usage: y = shoot(@f, ya, yb, z1, z2, t, tol)
  %
  % Returns:
  %     y:      array of solution function values corresponding to the
  %             values in the supplied array t.  The array y will be Nx1
  %             where N is the number of elements in the array t.
  %
  % Parameters:
  %     f:      handle of function dy/dt = f(y,t).  Since we are solving
  %             a second-order boundary-value problem that has been
  %             transformed into a first order system, this function
  %             should return a 1x2 array with the first entry equal
  %             to y and the second entry equal to y'.
  %     ya:     solution value at the left boundary: ya = y(a).
  %     yb:     solution value at the right boundary: yb = y(b).
  %     z1:     first estimate of y'(t(1)).
  %     z2:     second estimate of y'(t(1)).
  %     t:      array of time values from [a,b] at which to compute y.
  %     tol:    allowable tolerance on right boundary: |b - y(n)| < tol
  %
  % This function implements the shooting method to solve second-order boundary
  % value problems.  It assumes that the second order equation has been
  % converted to a first order system of two equations and uses the built-in
  % Octave function lsode() to solve the initial value problems.  The secant
  % method is used to refine the initial values of y' used for the initial
  % value problems.
  % --------------------------------------------------------------------------

  MAX_ITER = 25;   % Maximum number of shooting iterations

  n = max( size ( t ) ); % Determine the size of the arrays we will generate

  % Compute solution to first initial value problem (IVP) with y'(a) = z1.
  % Because we are using the secant method to refine our estimates of z = y',
  % we don't really need all the solution of the IVP, just the last point of
  % it -- this is saved in w1.

  y = lsode( f, [ya z1], t );
  w1 = y(n,1);

  % Begin the main loop.  We will compute the solution of a second IVP
  % and then use the both solutions to refine our estimate of y'(a).  This
  % second solution then replaces the first and a new "second" solution is
  % generated.  This process continues until we either solve the problem to
  % within the specified tolerance or we exceed the maximum number of
  % allowable iterations.

  for i = 1 : MAX_ITER

    % Solve second initial value problem, using y'(a) = z2.  We need to
    % retain the entire solution vector y since if y(t(n)) is close enough
    % to yb for us to stop then the first column of y becomes our solution
    % vector.

    y = lsode( f, [ya z2], t );
    w2 = y(n,1);

    fprintf( '%2d: z1 = %f, z2 = %f, error = %e\n', i, z1, z2, yb - w2 );
    fflush( stdout );

    % Check to see if we are done...
    
    if ( abs( yb - w2 ) < tol )
      break;
    end

    % Compute the new approximations to the initial value of the first
    % derivative.  We compute z2 using a linear fit through (z1,w1) and
    % (z2,w2) where w1 and w2 are the estimates at t=b of the initial
    % value problems solved above with y1'(a) = z1 and y2'(a) = z2.  The
    % new value for z1 is the old value of z2.  This is identical to the
    % secant method to solve the root-finding problem.

    temp = z2;
    z2 = z2 + ( z2 - z1 ) / ( w2 - w1 ) * ( yb - w2 );
    z1 = temp;
    w1 = w2;

  end

  % All done.  Check to see if we really solved the problem, and then return
  % the solution.

  if ( i >= MAX_ITER && abs( yb - w2 ) >= tol )
    fprintf( '\a**** ERROR ****\n' );
    fprintf( 'Maximum number of iterations (%d) exceeded\n', MAX_ITER );
    fprintf( 'Returned values may not have desired accuracy\n' );
    fprintf( 'Error estimate of returned solution is %e\n', yb - w2 );
  end

  y = y(:,1);

end

% End of file
