%
%   BVP2SH.M  This function file solves a 2nd order BVP with linear BCs
%                        via the Shooting Method
%
%   This method solves the BVP using an IVP solver with an initial guess
%   for the unknown initial conditions.  An iterative scheme automatically
%   modifies this guess until the boundary conditions at the second boundary
%   point are satisfied.  This iterative technique is called the Shooting Method.
%
%   A general 2nd order ODE can be written in the form
%        y'' + a*y' + b*y = f(t)
%   where, in general, a and b can be functions of y and t.
%
%   We can also write this as a set of two 1st order equations in so-called 
%   state space form as
%        dz1/dt = z2   and   dz2/dt = f - a*z2 - b*z1
%   where z1 = y  and  z2 = y'
%
%   Now, if we had a set of initial conditions, zo = [y(to); y'(to)], then
%   this state space system could be easily solved as an IVP.
%
%   However, since this is a BVP, we do not have a complete set of initial 
%   conditions.  Thus, the basic idea of the shooting method is to guess at the
%   unknown variable at the initial state so that the solution of the resulting 
%   IVP gives the desired final state for the BVP.  This is done iteratively by
%   converging on the free parameter, alf, that is included as part of the 
%   initial condition vector such that the desired final state is achieved.
%
%   For the case of a 2nd order 2-point BVP with linear BCs, a general way to  
%   write the boundary conditions is as follows:
%     at p1:  a1*y(p1) + b1*y'(p1) - w1 = 0
%     at p2:  a2*y(p2) + b2*y'(p2) - w2 = 0
%   The coefficients in these expressions are input to this function in the
%   zbc matrix as:  zbc = [a1 b1 w1; a2 b2 w2].  This format essentially 
%   defines the BCs for the 2nd order BVP with linear BCs (note that the equation 
%   can be nonlinear -- but not the BCs).
%
%   At the initial boundary point, p1, a check is made to determine the 
%   kind of boundary condition and thus the form of the initial condition 
%   vector for the iterative IVP solution:
%     Dirichlet BC:     if b1 = 0,  y(p1) = w1/a1   
%                       then      zo = [w1/a1; alf]          (guess initial slope)
%     Neumann BC:       if a1 = 0,  y'(p1) = w1/b1
%                       then      zo = [alf; w1/b1]          (guess initial value)
%     Robin (mixed) BC: if both a1 and b1 are nonzero, 
%                       then  zo = [alf; (w1 - a1*alf)/b1]   (guess initial value)
%   where alf is the unknown parameter which will be computed as part 
%   of the iterative solution.  alf is determined such that the second 
%   BC is satisfied to within some small tolerance. Note that the Neumann 
%   BC is just a subset of the Robin BC with a1 = 0, and it is implemented as
%   such in this routine.
%
%   The inputs to this function are similar to those used with Matlab's ode23
%   IVP solver (in its simplest implementation) except that the zbc matrix 
%   replaces the initial state vector.  With zbc known, the zo vector is computed 
%   internally as described above.  Also, there is an option to input a guess
%   for alf, if desired
%     func - name of function file that evaluates the RHS of the state equations
%     tspan   - vector that defines the range for the independent variable
%     zbc  - matrix containing the BC coefficients (as described above)
%     options - standard options structure for use in ode23  (optional)
%     alfo - inital guess for alf  (optional)
%
%   The outputs are:
%      t - column vector containing the independent variable 
%      z - solution matrix with two columns with column 1 containing z1(t)
%          and column 2 containing z2(t)
%
%   Note #1: Extra parameters that may be needed in the ode function file must be 
%   passed to the function file via a global statement in the main program and  
%   the function file.  The current, relatively simple implementation of the  
%   Shooting Method does not use the variable input arguments (varargin) that are  
%   often used in the ode23 routine.  This is not overly restrictive, however, since   
%   it is just as easy to pass needed variables into the function via the global 
%   statement.
%
%   Note #2: This routine assumes that you are integrating from boundary point
%   p1 to boundary point p2.  This normally implies that p2 > p1, but this is not 
%   absolutely essential.  In some cases, it may be appropriate to integrate 
%   'backwards'  -- just be careful with the ordering of the BCs if you try this...
%
%   Note #3:  Nonlinear problems are sometimes sensitive to the initial guess for
%   the free variable, alf, in the initial condition vector.  Here we have simply 
%   set alfo = 1 as default.  The user can override this as needed in particularly  
%   difficult problems.  A poor initial guess usually just means that a few more  
%   iterations may be needed. If the system does not converge within 10 - 20 
%   iterations, you may want to change alfo via the last entry in the input argument
%   list to see if this helps.  This is not usually necessary,  however...
%
%   File prepared by J. R. White, UMass-Lowell (March 2003)
%

      function [t,z]  = bvp2sh(func,tspan,zbc,options,alfo)
%
%   check to see if the options structure is available
      if nargin < 4,  options = [];   end
%
%   check to see if there is an initial guess for alfo
      if nargin < 5,  alfo = 1;   end
%
%   extract BC coefficients (for convenience)
      a1 = zbc(1,1);  b1 = zbc(1,2);  w1 = zbc(1,3);
      a2 = zbc(2,1);  b2 = zbc(2,2);  w2 = zbc(2,3);
%
%   set some other iterative parameters
      err = 1e10;   tol = 1e-6;   icnt = 1;   icntmax = 25;
%
%   start iteration  
      while abs(err) > tol  &  icnt <= icntmax
%
%   determine initial condition vector based on structure of zbc matrix
        if b1 == 0
          zo = [w1/a1; alfo];
        else
          zo = [alfo; (w1 - a1*alfo)/b1];
        end
%
%   solve IVP with current guess for initial conditions
        [t,z] = ode23(func,tspan,zo,options);
%
%   evaluate error in second BC
        e1 =  a2*z(end,1) + b2*z(end,2) - w2;  err = e1;
        fprintf(1,'  \n')
        fprintf(1,'For iteration # %2i: \n',icnt)
        fprintf(1,'  Initial conditions are: %13.5e %13.5e \n',zo)
        fprintf(1,'  Error in 2nd BC is:     %13.5e \n',err)
        fprintf(1,'  \n')
%
%   estimate new values for alf (if necessary)
        if abs(err) > tol
          alfp = 1.01*alfo;
          if b1 == 0
            zo = [w1/a1; alfp];
          else
            zo = [alfp; (w1 - a1*alfp)/b1];
          end
          [t,z] = ode23(func,tspan,zo,options);
          e2 =  a2*z(end,1) + b2*z(end,2) - w2;  
          deda = (e2-e1)/(0.01*alfo);    alfn = alfo-e1/deda;
          icnt = icnt+1;   alfo = alfn;
        end
      end   % end of iteration loop
      if icnt >= icntmax
        fprintf(1,'  WARNING -- Hit maximum iteration limit!!! \n');
      end
%
%   end of function