function ex8bvp
%EX8BVP  Example 8 of the BVP tutorial.
%   These are equations describing fluid injection through one side of
%   a long  vertical channel, considered in Example 1.4 of U.M. Ascher, 
%   R.M.M. Mattheij, and R.D. Russell, Numerical Solution of Boundary 
%   Value Problems for Ordinary Differential Equations, SIAM, 1995
%
%   The differential equations
%
%      f''' - R*[(f')^2 - f*f''] + R*A = 0
%      h'' + R*f*h' + 1 = 0
%      theta'' + P*f*theta' = 0
%
%   are to be solved subject to boundary conditions
%
%      f(0) = f'(0) = 0
%      f(1) = 1
%      f'(1) = 0
%      h(0) = h(1) = 0
%      theta(0) = 0
%      theta(1) = 1
%
%   Here R and P are known constants, but A is determined by the boundary 
%   conditions.
%
%   For a Reynolds number R = 100, this problem can be solved with crude
%   guesses, but as R increases, it becomes much more difficult because of
%   a boundary layer at x = 0.  The solution is computed here for R = 10000
%   by continuation, i.e., the solution for one value of R is used as guess
%   for R = R*10.  This is a comparatively expensive problem for BVP4C because
%   a fine mesh is needed to resolve the boundary layer and there are 7 
%   unknown functions and 1 unknown parameter.


  % The solution is first sought for R = 100.
  R = 100;

  % A crude mesh of 10 equally spaced points and a constant 1 
  % for each solution component is used for the first value of R.
  % The unknown parameter A is guessed to be 1.
  solinit = bvpinit(linspace(0,1,10),ones(7,1),1);

  sol = bvp4c(@ex8ode,@ex8bc,solinit);

  fprintf('For R = %5i, A = %4.2f.\n',R,sol.parameters);
  figure
  lines = {'k-.','r--','b-'};
  plot(sol.x,sol.y(2,:),lines{1});
  axis([-0.1 1.1 0 1.7]);
  title('Fluid injection problem');
  xlabel('x');
  ylabel('f'' (x)');
  drawnow
  
  % The solution is computed for larger R by continuation, i.e.,
  % the solution for one value of R is used as guess for the next.
  hold on
  for i=2:3
    R = R*10; 
    % After creating function handles, the new value of R 
    % will be used in nested functions.
    sol = bvp4c(@ex8ode,@ex8bc,sol);
    fprintf('For R = %5i, A = %4.2f.\n',R,sol.parameters);
    plot(sol.x,sol.y(2,:),lines{i});
    drawnow
  end
  legend('R =    100','R =   1000','R = 10000',1); 
  hold off

  % --------------------------------------------------------------------------
  % Nested functions
  %

  function dydx = ex8ode(x,y,A);
  %EX8ODE  ODE function for Example 8 of the BVP tutorial.
    P = 0.7*R;
    dydx = [ y(2)
             y(3)
             R*(y(2)^2 - y(1)*y(3) - A)
             y(5)
             -R*y(1)*y(5) - 1
             y(7)
             -P*y(1)*y(7) ];
  end % ex8ode
    
  % --------------------------------------------------------------------------

  function res = ex8bc(ya,yb,A)
  %EX8BC  Boundary conditions for Example 8 of the BVP tutorial.
    res = [ya(1)
           ya(2)
           yb(1) - 1
           yb(2)
           ya(4)
           yb(4)
           ya(6)
           yb(6) - 1];
  end % ex8bc

  % --------------------------------------------------------------------------
  
end  % ex8bvp