function [x,y,r] = bvp (alpha,beta,a,b,tol,m,q,g)
%----------------------------------------------------------------
% Usage:       [x,y,r] = bvp (alpha,beta,a,b,tol,m,q,g)
%
% Description: Use the finite difference method to solve the 
%              following boundary value problem where y' = dy/dx:
%
%              y'' = g(x,y,y')  ,  y(alpha) = a, y(beta) = b
%
% Inputs:      alpha = initial value of x
%              beta  = final value of x
%              a     = initial value of y,
%              b     = final value of y
%              tol   = upper bound on difference between estimates 
%              m     = number of interior solution points (m >= 1)
%              q     = maximum number of iterations (q >= 1)
%              g     = string containing name of a user-supplied 
%                      function which defines the second-order 
%                      ordinary differential equation to be solved.
%                      The form of g is:
%
%                      y2 = function g (x,y,y1)    
%
%                      When g is called it should return the value
%                      y2 = d2y/dx2 where y1 = dy/dx.
%                        
% Outputs:      x = (m+2) by 1 vector containing independent 
%                   variables uniformly distributed over [alpha,beta]
%               y = (m+2) by 1 vector containing solution points
%               r = number of iterations performed. If r < q, then
%                   the following termination criterion was 
%                   satisfied where Dy is the change in y between
%                   successive iterations:
%
%                   ||Dy|| < tol
%----------------------------------------------------------------

% Initialize

   tol = args (tol,0,tol,5,'bvp');
   m   = args (m,1,m,6,'bvp');
   q   = args (q,1,q,7,'bvp');

   k = 0;
   dy = 0.001;
   h = (beta - alpha)/(m + 1);
   y0 = zeros (m+2,1);
   z0 = zeros (m,1); 
   z1 = zeros (m,1); 
   J  = zeros (m,m);   
   T  = zeros (3,m);
   x = linspace (alpha,beta,m+2)';

% Solve differential equation 

   y = zeros (m+2,1);
   y(1) = a;
   y(m+2) = b;
   e = tol + 1;
   hwbar = waitbar(0,'Solving Boundary Value Problem: bvp');

   while (e >= tol) & (k < q)

% Evaluate Jacobian matrix using central differences 

      waitbar (max(k/q,tol/e))
      y0 = y;
      J = -2*eye(m,m);
      for j = 1 : m
         y0(j+1) = y(j+1) + dy;
         z0 = funbvp (y0,alpha,beta,m,g);
         y0(j+1) = y(j+1) - 2*dy;
         z1 = funbvp (y0,alpha,beta,m,g);
         y0(j+1) = y(j+1) + dy;
         for i = 1 : m
            if (j == i-1) | (j == i+1)
               J(i,j) = 1;
            end    
            J(i,j) = J(i,j) - h^2*(z1(i) - z0(i))/(2*dy);
         end
      end
            
% Solve the tridiagonal system: J*z1 = -z0 

      z0 = funbvp (y,alpha,beta,m,g);
      for i = 1 : m-1
         T(1,i) = J(i,i+1);
         T(2,i) = J(i,i);
         T(3,i) = J(i+1,i);
      end
      T(2,m) = J(m,m);
      [z1,Delta] = tridec (T,-z0);
      e = norm(z1,inf);
      for i = 1 : m
         y(i+1) = y(i+1) + z1(i);
      end 
      k = k + 1;
   end

   close(hwbar)
   r = k;   
   
function e = funbvp (z,alpha,beta,m,g)
%----------------------------------------------------------------
%
% Usage:       e = funbvp (z,alpha,beta,m,g);
%
% Description: This function is used by bvp.  Newton's vector 
%              method is used to find the (m=2) by 1 vector z 
%              such that the m by 1 vector e = 0.
%----------------------------------------------------------------

   chkvec (z,1,'funbvp');
   e = zeros (m,1);
   h = (beta - alpha)/(m + 1);
   for i = 2 : m+1
      x = alpha + (i-1)*h;
      dz = (z(i+1) - z(i-1))/(2*h);
      e(i-1) = z(i+1) - 2*z(i) + z(i-1) - h^2*feval(g,x,z(i),dz);
   end