用Matlab bvp4c 解带未知参数的常微分方程边值问题，怎么给出合适的初值？

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.
<p>% The solution is first sought for R = 100.<br />
R = 100;</p>
<p>options = [];    % place holder</p>
<p>% A crude mesh of 10 equally spaced points and a constant 1<br />
% for each solution component is used for the first value of R.<br />
% The unknown parameter A is guessed to be 1.<br />
solinit = bvpinit(linspace(0,1,10),ones(7,1),1);</p>
<p>sol = bvp4c(@ex8ode,@ex8bc,solinit,options,R);</p>
<p>fprintf('For R = %5i, A = %4.2f.\n',R,sol.parameters);<br />
figure<br />
lines = {'k-.','r--','b-'};<br />
plot(sol.x,sol.y(2,:),lines{1});<br />
axis([-0.1 1.1 0 1.7]);<br />
title('Fluid injection problem');<br />
xlabel('x');<br />
ylabel('f'' (x)');<br />
drawnow</p>
<p>% The solution is computed for larger R by continuation, i.e.,<br />
% the solution for one value of R is used as guess for the next.<br />
hold on<br />
for i=2:3<br />
R = R*10;<br />
sol = bvp4c(@ex8ode,@ex8bc,sol,options,R);<br />
fprintf('For R = %5i, A = %4.2f.\n',R,sol.parameters);<br />
plot(sol.x,sol.y(2,:),lines{i});<br />
drawnow<br />
end<br />
legend('R =    100','R =   1000','R = 10000',1);<br />
hold off</p>
<p>% --------------------------------------------------------------------------</p>
<p>function dydx = ex8ode(x,y,A,R);<br />
%EX8ODE  ODE function for Example 8 of the BVP tutorial.<br />
P = 0.7<em>R;<br />
dydx = [ y(2)<br />
y(3)<br />
R</em>(y(2)^2 - y(1)<em>y(3) - A)<br />
y(5)<br />
-R</em>y(1)<em>y(5) - 1<br />
y(7)<br />
-P</em>y(1)*y(7) ];</p>
<p>% --------------------------------------------------------------------------</p>
<p>function res = ex8bc(ya,yb,A,R)<br />
%EX8BC  Boundary conditions for Example 8 of the BVP tutorial.<br />
res = [ya(1)<br />
ya(2)<br />
yb(1) - 1<br />
yb(2)<br />
ya(4)<br />
yb(4)<br />
ya(6)<br />
yb(6) - 1];</p>

function ch3ex6
% CH3EX6 example 3.5.6 of book Solving ODEs with Matlab
<p>global R<br />
color = [ 'k', 'r', 'b' ];<br />
options = bvpset('FJacobian',@Jac,'BCJacobian',@BCJac,...<br />
'Vectorized','on');<br />
R = 100;<br />
sol = bvpinit(linspace(0,1,10),ones(7,1),1);<br />
hold on<br />
for i = 1:3<br />
sol = bvp4c(@ode,@bc,sol,options);<br />
fprintf('For R = %5i, A = %4.2f.\n',R,sol.parameters);<br />
plot(sol.x,sol.y(2,:),color(i));<br />
axis([-0.1 1.1 0 1.7]);<br />
drawnow<br />
R = 10<em>R;<br />
end<br />
legend('R = 100','R = 1000','R = 10000');<br />
hold off<br />
%================================================================<br />
function dydx = ode(x,y,A);<br />
global R<br />
P = 0.7</em>R;<br />
dydx = [ y(2,:); y(3,:); R<em>(y(2,:).^2- y(1,:).</em>y(3,:)-A);...<br />
y(5,:); -R<em>y(1,:).</em>y(5,:)-1; y(7,:); -P<em>y(1,:).</em>y(7,:) ];<br />
function [dFdy,dFdA] = Jac(x,y,A)<br />
global R</p>
<p>dFdy = [ 0, 1, 0, 0, 0, 0, 0<br />
0, 0, 1, 0, 0, 0, 0<br />
-R<em>y(3), 2</em>R<em>y(2), -R</em>y(1), 0, 0, 0, 0<br />
0, 0, 0, 0, 1, 0, 0<br />
-R<em>y(5), 0, 0, 0, -R</em>y(1), 0, 0<br />
0, 0, 0, 0, 0, 0, 1<br />
-7/10<em>R</em>y(7), 0, 0, 0, 0, 0, -7/10<em>R</em>y(1) ];<br />
dFdA=[0;0;-R;0;0;0;0];<br />
function res = bc(ya,yb,A)<br />
res = [ ya(1); ya(2); yb(1)-1; yb(2);...<br />
ya(4); yb(4); ya(6); yb(6)-1 ];<br />
function [dBCdya,dBCdyb,dBCdA] = BCJac(ya,yb,A)<br />
dBCdya=[1,0,0,0,0,0,0<br />
0, 1, 0, 0, 0, 0, 0<br />
0, 0, 0, 0, 0, 0, 0<br />
0, 0, 0, 0, 0, 0, 0<br />
0, 0, 0, 1, 0, 0, 0<br />
0, 0, 0, 0, 0, 0, 0<br />
0, 0, 0, 0, 0, 1, 0<br />
0, 0, 0, 0, 0, 0, 0 ];<br />
dBCdyb=[0,0,0,0,0,0,0<br />
0, 0, 0, 0, 0, 0, 0<br />
1, 0, 0, 0, 0, 0, 0<br />
0, 1, 0, 0, 0, 0, 0<br />
0, 0, 0, 0, 0, 0, 0<br />
0, 0, 0, 1, 0, 0, 0<br />
0, 0, 0, 0, 0, 0, 0<br />
0, 0, 0, 0, 0, 1, 0 ];<br />
dBCdA = zeros(8,1);</p>