Bvp4c_Exmamples.

txt
function bratubvp
%BRATUBVP Exercise for Example 1 of the BVP tutorial.
% The BVP y'' + exp(y) = 0, y(0) = 0 = y(1) is a
standard example
% of a problem with two solutions. It is easy enough to
solve, but
% some experimentation with the guess may be necessary to
get both.

% Copyright 1999, The MathWorks, Inc.

options = bvpset('stats','on');
solinit = bvpinit(linspace(0,1,5),[0.1 0]);
sol1 = bvp4c(@bratuode,@bratubc,solinit,options);
fprintf('\n');
% Change the initial guess to converge to a different
solution.
solinit = bvpinit(linspace(0,1,5),[3 0]);
sol2 = bvp4c(@bratuode,@bratubc,solinit,options);

clf reset
plot(sol1.x,sol1.y(1,:),sol2.x,sol2.y(1,:))
title('Bratu''s equation has two solutions when \lambda =
1.')
xlabel('x')
ylabel('y')
shg
%
-----------------------------------------------------------
---------------

function dydx = bratuode(x,y)

%BRATUODE ODE function for the exercise of Example 1 of
the BVP tutorial.
dydx = [ y(2)
-exp(y(1))];

%
-----------------------------------------------------------
---------------

Page 1
 Bvp4c_Exmamples.txt
function res = bratubc(ya,yb)
%BRATUBC Boundary conditions for the exercise of Example 1
of the BVP tutorial.
res = [ya(1)
yb(1)];

function ex1bvp
%EX1BVP Example 1 of the BVP tutorial.
% This is the example for MUSN in U. Ascher, R. Mattheij,
and R. Russell,
% Numerical Solution of Boundary Value Problems for
Ordinary Differential
% Equations, SIAM, Philadelphia, PA, 1995. MUSN is a
multiple shooting
% code for nonlinear BVPs. The problem is
%
% u' = 0.5*u*(w - u)/v
% v' = -0.5*(w - u)
% w' = (0.9 - 1000*(w - y) - 0.5*w*(w - u))/z
% z' = 0.5*(w - u)
% y' = -100*(y - w)
%
% The interval is [0 1] and the boundary conditions are
%
% u(0) = v(0) = w(0) = 1, z(0) = -10, w(1) = y(1)
%
% The example uses a guess for the solution coded here in
EX1INIT.
% The results of a run of the FORTRAN code MUSN are here
compared to
% the curves produced by BVP4C.

% Copyright 2002, The MathWorks, Inc.

solinit = bvpinit(linspace(0,1,5),@ex1init);
options = bvpset('Stats','on','RelTol',1e-5);
sol = bvp4c(@ex1ode,@ex1bc,solinit,options);

% The solution at the mesh points

x = sol.x;
y = sol.y;

% Solution obtained using MUSN:

Page 2
 Bvp4c_Exmamples.txt
amrx = [ 0. .1 .2 .3 .4 .5 .6 .7 .8 .9 1.]';
amry = [1.00000e+00 1.00000e+00 1.00000e+00
-1.00000e+01 9.67963e-01
1.00701e+00 9.93036e-01 1.27014e+00
-9.99304e+00 1.24622e+00
1.02560e+00 9.75042e-01 1.47051e+00
-9.97504e+00 1.45280e+00
1.05313e+00 9.49550e-01 1.61931e+00
-9.94955e+00 1.60610e+00
1.08796e+00 9.19155e-01 1.73140e+00
-9.91915e+00 1.72137e+00
1.12900e+00 8.85737e-01 1.81775e+00
-9.88574e+00 1.80994e+00
1.17554e+00 8.50676e-01 1.88576e+00
-9.85068e+00 1.87957e+00
1.22696e+00 8.15025e-01 1.93990e+00
-9.81503e+00 1.93498e+00
1.28262e+00 7.79653e-01 1.98190e+00
-9.77965e+00 1.97819e+00
1.34161e+00 7.45374e-01 2.01050e+00
-9.74537e+00 2.00827e+00
1.40232e+00 7.13102e-01 2.02032e+00
-9.71310e+00 2.02032e+00];
% Shift up the fourth component for the plot.
amry(:,4) = amry(:,4) + 10;
y(4,:) = y(4,:) + 10;
figure
plot(x,y',amrx,amry,'*')
axis([0 1 -0.5 2.5])
title('Example problem for MUSN')
ylabel('bvp4c and MUSN (*) solutions')
xlabel('x')

%
-----------------------------------------------------------
---------------

function dydx = ex1ode(x,y)

%EX1ODE ODE function for Example 1 of the BVP tutorial.
% The components of y correspond to the original
variables
% as y(1) = u, y(2) = v, y(3) = w, y(4) = z, y(5) = y.
Page 3
 Bvp4c_Exmamples.txt
dydx = [ 0.5*y(1)*(y(3) - y(1))/y(2)
-0.5*(y(3) - y(1))
(0.9 - 1000*(y(3) - y(5)) - 0.5*y(3)*(y(3) -
y(1)))/y(4)
0.5*(y(3) - y(1))
100*(y(3) - y(5)) ];

%----------------------------------------------------------
---------------
function res = ex1bc(ya,yb)
%EX1BC Boundary conditions for Example 1 of the BVP
tutorial.
% RES = EX1BC(YA,YB) returns a column vector RES of the
% residual in the boundary conditions resulting from the
% approximations YA and YB to the solution at the ends of
% the interval [a b]. The BVP is solved when RES = 0.
% The components of y correspond to the original
variables
% as y(1) = u, y(2) = v, y(3) = w, y(4) = z, y(5) = y.
res = [ ya(1) - 1
ya(2) - 1
ya(3) - 1
ya(4) + 10
yb(3) - yb(5)];
%----------------------------------------------------------
---------------
function v = ex1init(x)
%EX1INIT Guess for the solution of Example 1 of the BVP
tutorial.
v = [ 1
1
-4.5*x^2+8.91*x+1
-10
-4.5*x^2+9*x+0.91 ];

Page 4
 Bvp4c_Exmamples.txt

function ex2bvp
%EX2BVP Example 2 of the BVP tutorial.
% A standard linear problem with a boundary layer at the
origin.
% The differential equation y'' + 3*p*y/(p + t^2)^2 = 0
has the
% analytical solution y(t) = t/sqrt(p + t^2). The
parameter p
% is taken to be 1e-5, a common value in tests. The
solution is
% to have specified values at t = -0.1 and +0.1, values
taken from
% this analytical solution.
%
% The default RelTol of 1e-3 gives an acceptable
solution, but
% reducing RelTol to 1e-4 resolves better the boundary
layer. A
% constant guess is used for RelTol = 1e-3. The same
guess could be
% used for RelTol = 1e-4, but a very much better guess is
provided
% by the solution previously computed for RelTol = 1e-3.
% Copyright 2002, The MathWorks, Inc.

% Evaluate the analytical solution for comparison.

tt = -0.1:0.01:+0.1;
p = 1e-5;
yy = tt ./ sqrt(p + tt .^2);
options = bvpset('stats','on','Fjacobian',@ex2Jac);

% BVPINT is used to specify an initial guess for the mesh

of 10
% equally spaced points. A constant guess based on a
straight line
% between the boundary values for y is 0 for y(t) and 10
Page 5
 Bvp4c_Exmamples.txt
for y'(t).
solinit = bvpinit(linspace(-0.1,0.1,10),[0 10]);
sol = bvp4c(@ex2ode,@ex2bc,solinit, options);
t = sol.x;
y = sol.y;

figure
plot(t,y(1,:),tt,yy,'*')
axis([-0.1 0.1 -1.1 1.1])
title(['Linear boundary layer problem with RelTol =
1e-3.'])
xlabel('t')
ylabel('y and analytical (*) solutions')
fprintf('\n');
% A smaller RelTol is used to resolve better the boundary
layer.
% The previous solution provides an excellent guess.
options = bvpset(options,'RelTol',1e-4);
sol = bvp4c(@ex2ode,@ex2bc,sol,options);
t = sol.x;
y = sol.y;

figure
plot(t,y(1,:),tt,yy,'*')
axis([-0.1 0.1 -1.1 1.1])
title(['Linear boundary layer problem with RelTol =
1e-4.'])
xlabel('t')
ylabel('y and analytical (*) solutions')

%
-----------------------------------------------------------
---------------
function dydt = ex2ode(t,y)
%EX2ODE ODE function for Example 2 of the BVP tutorial.
% The components of y correspond to the original
variables
% as y(1) = y, y(2) = y'.
p = 1e-5;
dydt = [ y(2)
Page 6
 Bvp4c_Exmamples.txt
-3*p*y(1)/(p+t^2)^2];
%
-----------------------------------------------------------
---------------

function dfdy = ex2Jac(t,y)

%EX2JAC The Jacobian of the ODE function for Example 2 of
the BVP tutorial.
p = 1e-5;
dfdy = [ 0 1
-3*p/(p+t^2)^2 0 ];

%
-----------------------------------------------------------
---------------
function res = ex2bc(ya,yb)
%EX2BC Boundary conditions for Example 2 of the BVP
tutorial.
% The boundary conditions are that the solution should
agree
% with the values of an analytical solution at both a and
b.
p = 1e-5;
yatb = 0.1/sqrt(p + 0.01);
yata = - yatb;
res = [ ya(1) - yata
yb(1) - yatb ]; function ex3bvp
%EX3BVP Example 3 of the BVP tutorial.
% This is the example for D02KAF from the NAG library.
D02KAF is a code
% for Sturm-Liouville problems that takes advantage of
special features of
% the problem. The task is to find the fourth
eigenvalue of Mathieu's
% equation
%
% y'' + (lambda -2*q*cos(2*x))*y = 0
%
% on the interval [0, pi] with boundary conditions y'(0)
= 0, y'(pi) = 0.
% The parameter q = 5.
%
Page 7
 Bvp4c_Exmamples.txt
% A code that exploits fully the special nature of the
Sturm-Liouville
% problem can compute a specific eigenvalue. Of course
using BVP4C we can
% only compute an eigenvalue near to a guessed value. We
can make it much
% more likely that we compute the desired eigenfunction
by supplying a
% guess that has the correct qualitative behavior. We
use here the same
% guess lambda = 15 as the example. The eigenfunction is
determined only
% to a constant multiple, so the normalizing condition
y(0) = 1 is used to
% specify a particular solution.
%
% Plotting the solution on the mesh found by BVP4C does
not result in a
% smooth graph. The solution S(x) is continuous and has
a continuous
% derivative. It can be evaluated inexpensively using
DEVAL at as many
% points as necessary to get a smooth graph.
% Copyright 2004, The MathWorks, Inc.

% BVPINT is used to form an initial guess for a mesh of 10

equally
% spaced points. The guess cos(4x) for y(x) and its
derivative as guess
% for y'(x) are evaluated in EX3INIT. The desired
eigenvalue is the one
% nearest the guess lambda = 15. A guess for unknown
parameters is the
% (optional) last argument of BVPINIT.
solinit = bvpinit(linspace(0,pi,10),@ex3init,15);
options = bvpset('stats','on');

sol = bvp4c(@ex3ode,@ex3bc,solinit,options);
% BVP4C returns the solution as the structure 'sol'. The
computed eigenvalue
% is returned in the field sol.parameters.
Page 8
 Bvp4c_Exmamples.txt
fprintf('\n');
fprintf('D02KAF computed lambda = 17.097.\n')
fprintf('bvp4c computed lambda =%7.3f.\n',sol.parameters)

figure
plot(sol.x,sol.y(1,:),sol.x,sol.y(1,:),'*')
axis([0 pi -1 1])
title('Eigenfunction for Mathieu''s equation.')
xlabel('Solution at mesh points only.')
% Plotting the solution just at the mesh points does not
result in a
% smooth graph near the ends of the interval. The
approximate solution
% S(x) is continuous and has a continuous derivative.
DEVAL is used to
% evaluate it at enough points to get a smooth graph.
figure
xint = linspace(0,pi);
Sxint = deval(sol,xint);
plot(xint,Sxint(1,:))
axis([0 pi -1 1])
title('Eigenfunction for Mathieu''s equation.')
xlabel('Solution evaluated on a finer mesh with DEVAL.')

%
-----------------------------------------------------------
---------------

function dydx = ex3ode(x,y,lambda)

%EX3ODE ODE function for Example 3 of the BVP tutorial.
q = 5;
dydx = [ y(2)
-(lambda - 2*q*cos(2*x))*y(1) ];

%
-----------------------------------------------------------
---------------

function res = ex3bc(ya,yb,lambda)

%EX3BC Boundary conditions for Example 3 of the BVP
tutorial.
res = [ ya(2)
yb(2)
Page 9
 Bvp4c_Exmamples.txt
ya(1) - 1 ];
%
-----------------------------------------------------------
---------------

function v = ex3init(x)
%EX3INIT Guess for the solution of Example 3 of the BVP
tutorial.
v = [ cos(4*x)
-4*sin(4*x) ];

function ex3bvp
%EX3BVP Example 3 of the BVP tutorial.
% This is the example for D02KAF from the NAG library.
D02KAF is a code
% for Sturm-Liouville problems that takes advantage of
special features of
% the problem. The task is to find the fourth
eigenvalue of Mathieu's
% equation
%
% y'' + (lambda -2*q*cos(2*x))*y = 0
%
% on the interval [0, pi] with boundary conditions y'(0)
= 0, y'(pi) = 0.
% The parameter q = 5.
%
% A code that exploits fully the special nature of the
Sturm-Liouville
% problem can compute a specific eigenvalue. Of course
using BVP4C we can
% only compute an eigenvalue near to a guessed value. We
can make it much
% more likely that we compute the desired eigenfunction
by supplying a
% guess that has the correct qualitative behavior. We
use here the same
% guess lambda = 15 as the example. The eigenfunction is
determined only
% to a constant multiple, so the normalizing condition
y(0) = 1 is used to
Page 10
 Bvp4c_Exmamples.txt
% specify a particular solution.
%
% Plotting the solution on the mesh found by BVP4C does
not result in a
% smooth graph. The solution S(x) is continuous and has
a continuous
% derivative. It can be evaluated inexpensively using
DEVAL at as many
% points as necessary to get a smooth graph.
% Copyright 2004, The MathWorks, Inc.

% BVPINT is used to form an initial guess for a mesh of 10

equally
% spaced points. The guess cos(4x) for y(x) and its
derivative as guess
% for y'(x) are evaluated in EX3INIT. The desired
eigenvalue is the one
% nearest the guess lambda = 15. A guess for unknown
parameters is the
% (optional) last argument of BVPINIT.
solinit = bvpinit(linspace(0,pi,10),@ex3init,15);
options = bvpset('stats','on');

sol = bvp4c(@ex3ode,@ex3bc,solinit,options);
% BVP4C returns the solution as the structure 'sol'. The
computed eigenvalue
% is returned in the field sol.parameters.
fprintf('\n');
fprintf('D02KAF computed lambda = 17.097.\n')
fprintf('bvp4c computed lambda =%7.3f.\n',sol.parameters)

figure
plot(sol.x,sol.y(1,:),sol.x,sol.y(1,:),'*')
axis([0 pi -1 1])
title('Eigenfunction for Mathieu''s equation.')
xlabel('Solution at mesh points only.')

% Plotting the solution just at the mesh points does not

result in a
% smooth graph near the ends of the interval. The
approximate solution
Page 11
 Bvp4c_Exmamples.txt
% S(x) is continuous and has a continuous derivative.
DEVAL is used to
% evaluate it at enough points to get a smooth graph.
figure
xint = linspace(0,pi);
Sxint = deval(sol,xint);
plot(xint,Sxint(1,:))
axis([0 pi -1 1])
title('Eigenfunction for Mathieu''s equation.')
xlabel('Solution evaluated on a finer mesh with DEVAL.')
%
-----------------------------------------------------------
---------------
function dydx = ex3ode(x,y,lambda)
%EX3ODE ODE function for Example 3 of the BVP tutorial.
q = 5;
dydx = [ y(2)
-(lambda - 2*q*cos(2*x))*y(1) ];
%
-----------------------------------------------------------
---------------

function res = ex3bc(ya,yb,lambda)

%EX3BC Boundary conditions for Example 3 of the BVP
tutorial.
res = [ ya(2)
yb(2)
ya(1) - 1 ];

%
-----------------------------------------------------------
---------------
function v = ex3init(x)
%EX3INIT Guess for the solution of Example 3 of the BVP
tutorial.
v = [ cos(4*x)
-4*sin(4*x) ];
function ex4bvp
%EX4BVP Example 4 of the BVP tutorial.
Page 12
 Bvp4c_Exmamples.txt
% This is a nerve impulse model considered in Example 7.1
of
% R. Seydel, From Equilibrium to Chaos, Elsevier, 1988.
The
% differential equations
%
% y1' = 3*(y1 + y2 - 1/3*y1^3 + lambda)
% y2' = -(y1 - 0.7 + 0.8*y2)/3
%
% are to be solved subject to periodic boundary
conditions
%
% y1(0) = y1(T)
% y2(0) = y2(T)
%
% This is an example of non-separated boundary
conditions, meaning
% that some conditions involve values of the solution at
both ends
% of the interval. Periodic boundary conditions are a
common source
% of such problems. BVP4C is unusual in that it accepts
problems
% with non-separated boundary conditions.
%
% The parameter lambda has the value -1.3 in the text
cited. The
% period T is unknown, which is to say that the length of
the
% interval is not known. Such problems require some
preparation.
% By scaling the independent variable to tau = t/T the
problem is
% posed on the fixed interval [0,1]. When this is done, T
becomes
% an unknown parameter, but BVP4C provides for unknown
parameters.

% Copyright 2002, The MathWorks, Inc.

% Periodic functions evaluated in EX4INIT are used as a

guess for the
% solution on a crude mesh of 5 equally spaced points. A
guess of 2*pi
Page 13
 Bvp4c_Exmamples.txt
% is provided for the unknown parameter T.
solinit = bvpinit(linspace(0,1,5),@ex4init,2*pi);
options = bvpset('stats','on');

sol = bvp4c(@ex4ode,@ex4bc,solinit,options);
T = sol.parameters;

% The independent variable needs to be rescaled to its

original value
% before plotting the solution. The initial guess is also
plotted.
figure
plot(T*sol.x,sol.y(1,:),T*solinit.x,solinit.y(1,:),'o')
legend('solution found','initial guess');
axis([0 T -2.2 2]);
title('Nerve impulse model');
ylabel('solution y_1(t)');
xlabel(['Periodic solution with period
',num2str(sol.parameters,4)]);
%
-----------------------------------------------------------
---------------

function v = ex4init(x)
%EX4INIT Guess function for Example 4 of the BVP tutorial.
% V = EX4INIT(X) returns a column vector V that is a
guess for y(x).
v = [ sin(2*pi*x)
cos(2*pi*x)];

%
-----------------------------------------------------------
---------------
function dydt = ex4ode(t,y,T);
%EX4ODE ODE funcion for Example 4 of the BVP tutorial.
dydt = [ 3*T*(y(1) + y(2) - 1/3*(y(1)^3) - 1.3);
(-1/3)*T*(y(1) - 0.7 + 0.8*y(2)) ];
%
-----------------------------------------------------------
---------------
Page 14
 Bvp4c_Exmamples.txt

function res = ex4bc(ya,yb,T)

%EX4BC Boundary conditions for Example 4 of the BVP
tutorial.
res = [ya(1) - yb(1)
ya(2) - yb(2)
T*(-1/3)*(ya(1) - 0.7 + 0.8*ya(2)) - 1];

function ex5bvp
%EX5BVP Example 5 of the BVP tutorial.
% Falkner-Skan BVPs are discussed in T. Cebeci and H.B.
Keller,
% Shooting and parallel shooting methods for solving the
Falkner-Skan
% boundary-layer equation, J. Comp. Phy., 7 (1971)
289-300. This is
% the positive wall shear case for which the parameter
beta is known
% and the problem is to be solved for a range of the
parameter. This
% is the hardest case of the table in the paper.
%
% The problem is posed on [0 infinity). As in the paper
cited, the
% boundary condition at infinity is imposed at a finite
point, here
% called 'infinity'. It is best to start with a
relatively small value
% and increase it until consistent results are obtained.
A value of 6
% appears to be satisfactory. Starting with a "large"
value is tempting,
% but not a good tactic because the code will fail with
the crude guess
% and default tolerances used here.
% Copyright 1999, The MathWorks, Inc.

% 'infinity' is a variable to facilitate experimentation.

infinity = 6;
% The constant guess for the solution satisfies the
boundary conditions.
solinit = bvpinit(linspace(0,infinity,5),[0 0 1]);
Page 15
 Bvp4c_Exmamples.txt

options = bvpset('stats','on');
sol = bvp4c(@ex5ode,@ex5bc,solinit,options);
eta = sol.x;
f = sol.y;

fprintf('\n');
fprintf('Cebeci & Keller report f''''(0) = 0.92768.\n')
fprintf('Value computed here is f''''(0) =
%7.5f.\n',f(3,1))

clf reset
plot(eta,f(2,:));
axis([0 infinity 0 1.4]);
title('Falkner-Skan equation, positive wall shear, \beta =
0.5.')
xlabel('\eta')
ylabel('df/d\eta')
shg
%
-----------------------------------------------------------
---------------

function dfdeta = ex5ode(eta,f)

%EX5ODE ODE function for Example 5 of the BVP tutorial.
beta = 0.5;
dfdeta = [ f(2)
f(3)
-f(1)*f(3) - beta*(1 - f(2)^2) ];

%
-----------------------------------------------------------
---------------
function res = ex5bc(f0,finf)
%EX5BC Boundary conditions for Example 5 of the BVP
tutorial.
res = [f0(1)
f0(2)
finf(2) - 1];

function ex6bvp
Page 16
 Bvp4c_Exmamples.txt
%EX6BVP Example 6 of the BVP tutorial.
% This is Example 2 from M. Kubicek et alia, Test
examples for comparison
% of codes for nonlinear boundary value problems in
ordinary differential
% equations, B. Childs et al., eds., Codes for
Boundary-Value Problems in
% Ordinary Differential Equations, Lecture Notes in
Computer Science #76,
% Springer, New York, 1979. This example shows how to
deal with a singular
% coefficient arising from reduction of a partial
differential equation to
% an ODE by symmetry. Also, for the physical parameters
considered here,
% the problem has three solutions.
% Copyright 2002, The MathWorks, Inc.

% Define the physical parameters for this problem.

f = 0.6;
g = 40;
b = 0.2;
options = []; % place holder

guess = [1; 0.5; 0];

others = [0.9070; 0.3639; 0.0001]; % Values reported
by Kubicek et alia.
fprintf('y(0): bvp4c Kubicek et al.\n')
figure

for index = 1:3

solinit = bvpinit(linspace(0,1,5),[guess(index) 0]);
sol = bvp4c(@ex6ode,@ex6bc,solinit,options,f,g,b);
fprintf(' %6.4f
%6.4f\n',sol.y(1,1),others(index))
plot(sol.x,sol.y(1,:))
if index == 1
axis([0 1 -0.1 1.1])
title('Multiple solutions to spherical catalyst
problem.')
Page 17
 Bvp4c_Exmamples.txt
xlabel('x')
ylabel('y')
hold on
end
drawnow
end
hold off

%
-----------------------------------------------------------
---------------

function dydx = ex6ode(x,y,f,g,b)

%EX6ODE ODE function for Example 6 of the BVP tutorial.
dydx = [y(2); 0];
temp = f^2 * y(1) * exp(g*b*(1-y(1))/(1+b*(1-y(1))));
if x == 0
dydx(2) = (1/3)*temp;
else
dydx(2) = -(2/x)*y(2) + temp;
end

%
-----------------------------------------------------------
---------------

function res = ex6bc(ya,yb,f,g,b)

%EX6BC Boundary conditions for Example 6 of the BVP
tutorial.
res = [ ya(2)
yb(1) - 1];

function ex7bvp
%EX7BVP Example 7 of the BVP tutorial.
% This is the first example problem for D02HBF from the
NAG library.
% The problem is y'' = (y^3 - y')/(2x), y(0) = 0.1, y(16)
= 1/6. The
% singularity at the origin is handled by using a series
to represent
% the solution and its derivative at a "small" distance d
> 0, namely
%
% y(d) = 0.1 + y'(0)*sqrt(d)/10 + d/100
Page 18
 Bvp4c_Exmamples.txt
% y'(d) = y'(0)/(20*sqrt(d)) + 1/100
%
% The value y'(0) is treated as an unknown parameter p.
The problem is
% solved numerically on [d, 16]. Two boundary conditions
are that the
% computed solution and its first derivative agree with
the values from
% the series at d. The remaining boundary condition is
y(16) = 1/6.
%
% The results from the documentation for D02HBF are here
compared to
% curves produced by bvp4c.
% Copyright 1999, The MathWorks, Inc.
options = bvpset('stats','on');

d = 0.1; % The known parameter

solinit = bvpinit(linspace(d,16,5),[1 1],0.2);

sol = bvp4c(@ex7ode,@ex7bc,solinit,options,d);
% Augment the solution array with the values y(0) = 0.1,
y'(0) = p
% to get a solution on [0, 16].
x = [0 sol.x];
y = [[0.1; sol.parameters] sol.y];
% Solution obtained using D02HBF, augmented with the values
y(0), y'(0).
nagx = [0.00 0.10 3.28 6.46 9.64 12.82 16.00]';
nagy = [0.1000 4.629e-2
0.1025 0.0173
0.1217 0.0042
0.1338 0.0036
0.1449 0.0034
0.1557 0.0034
0.1667 0.0035];
clf reset
plot(x,y(1,:),x,y(2,:),'--',nagx,nagy,'*')
axis([-1 16 0 0.18])
Page 19
 Bvp4c_Exmamples.txt
title('Problem with singular behavior at the origin.')
xlabel('x')
ylabel('y, dy/dx (--), and D02HBF (*) solutions')
shg

%
-----------------------------------------------------------
---------------

function dydx = ex7ode(x,y,p,d)

%EX7ODE ODE function for Example 7 of the BVP tutorial.
dydx = [ y(2)
(y(1)^3 - y(2))/(2*x) ];

%
-----------------------------------------------------------
---------------

function res = ex7bc(ya,yb,p,d)

%EX7BC Boundary conditions for Example 7 of the BVP
tutorial.
% The boundary conditions at x = d are that y and y' have
% values yatd and ypatd obtained from series expansions.
% The unknown parameter p = y'(0) and known parameter d
are
% used in the expansions.
yatd = 0.1 + p*sqrt(d)/10 + d/100;
ypatd = p/(20*sqrt(d)) + 1/100;
res = [ ya(1) - yatd
ya(2) - ypatd
yb(1) - 1/6 ];
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
%
Page 20
 Bvp4c_Exmamples.txt
% 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.
% Copyright 1999, The MathWorks, Inc.

% The solution is first sought for R = 100.

R = 100;
options = []; % place holder

% 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);
Page 21
 Bvp4c_Exmamples.txt

sol = bvp4c(@ex8ode,@ex8bc,solinit,options,R);
fprintf('For R = %5i, A = %4.2f.\n',R,sol.parameters);
clf reset
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;
sol = bvp4c(@ex8ode,@ex8bc,sol,options,R);
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
%
-----------------------------------------------------------
---------------

function dydx = ex8ode(x,y,A,R);

%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) ];

%
Page 22
 Bvp4c_Exmamples.txt
-----------------------------------------------------------
---------------
function res = ex8bc(ya,yb,A,R)
%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

function ex9bvp
%EX9BVP Example 9 of the BVP tutorial.
% This boundary value problem is the subject of Chapter 8
of
% C.C. Lin and L.A. Segel, Mathematics Applied to
Deterministic
% Problems in the Natural Sciences, SIAM, Philadelphia,
1988.
% The ODEs
%
% v' = (C - 1)/n
% C' = (vC - min(x,1))/eta
%
% are solved on the interval [0, lambda]. The boundary
conditions
% are v(0) = 0, C(lambda) = 1, and continuity of v(x) and
C(x) at
% x = 1. Accordingly, this is a three-point BVP that
must be
% reformulated for solution with the two-point BVP solver
BVP4C.
% This reformulation involves introducing unknowns y_1(x)
Page 23
 Bvp4c_Exmamples.txt
for v
% and y_2(x) for C on the interval 0 <= x <= 1 and
unknowns y_3(x)
% for v and y_4(x) for C on 1 <= x <= lambda. A new
independent
% variable is introduced for the second interval, tau =
% (x - 1)/(lambda - 1), so that it also ranges from 0 to
1. The
% differential equations for the four unknowns are then
solved
% on the interval [0, 1]. The continuity conditions on v
and C
% become boundary conditions on the new unknowns. A plot
of v(x)
% and C(x) on [0, lambda] involves plotting the new
unknowns over
% the subintervals.
%
% The quantity of most interest is the emergent
osmolarity Os =
% 1/v(lambda). The parameters are related to another
parameter
% kappa by eta = lambda^2/(n*kappa^2). Lin and Segel
develop an
% approximate solution for Os valid for "small" n. Here
the BVP is
% solved for a range of kappa when lambda = 2 and n =
0.005. The
% computed Os is compared to the approximation of Lin and
Segel.
% Copyright 2002, The MathWorks, Inc.

n = 5e-2;
lambda = 2;
options = []; % place holder

sol = bvpinit(linspace(0,1,5),[1 1 1 1]);

fprintf(' kappa computed Os approximate Os \n')

for kappa = 2:5
eta = lambda^2/(n*kappa^2);

Page 24
 Bvp4c_Exmamples.txt
sol = bvp4c(@ex9ode,@ex9bc,sol,options,n,lambda,eta);
K2 = lambda*sinh(kappa/lambda)/(kappa*cosh(kappa));
approx = 1/(1 - K2);
computed = 1/sol.y(3,end);
fprintf(' %2i %10.3f %10.3f
\n',kappa,computed,approx);
end

% v and C are computed separately on 0 <= x <= 1 and 1 <= x

<= lambda.
% A change of independent variable is used for the second
interval,
% which must then be undone to obtain the corresponding
mesh.
x = [sol.x sol.x*(lambda-1)+1];
y = [sol.y(1:2,:) sol.y(3:4,:)];

figure
plot(x,y(1,:),x,y(2,:),'--')
legend('v(x)','C(x)')
title('A three-point BVP.')
xlabel(['\lambda = ',num2str(lambda),', \kappa =
',num2str(kappa),'.'])
ylabel('v and C')

%
-----------------------------------------------------------
---------------
function dydx = ex9ode(x,y,n,lambda,eta)
%EX9ODE ODE function for Example 9 of the BVP tutorial.
dydx = [ (y(2) - 1)/n
(y(1)*y(2) - x)/eta
(lambda - 1)*(y(4) - 1)/n
(lambda - 1)*(y(3)*y(4) - 1)/eta ];
%
-----------------------------------------------------------
---------------
function res = ex9bc(ya,yb,n,lambda,eta)
%EX9BC Boundary conditions for Example 9 of the BVP
tutorial.
Page 25
 Bvp4c_Exmamples.txt
res = [ ya(1)
yb(4) - 1
yb(1) - ya(3)
yb(2) - ya(4)];

function ex9mbvp
%EX9MBVP Example 9 of the BVP tutorial, solved as a
multi-point BVP
% This boundary value problem is the subject of Chapter 8
of
% C.C. Lin and L.A. Segel, Mathematics Applied to
Deterministic
% Problems in the Natural Sciences, SIAM, Philadelphia,
1988.
% The ODEs
%
% v' = (C - 1)/n
% C' = (vC - min(x,1))/eta
%
% are solved on the interval [0, lambda]. The boundary
conditions
% are v(0) = 0, C(lambda) = 1, and continuity of v(x) and
C(x) at
% x = 1. Example EX9BVP shows how this three-point BVP
is
% reformulated for solution with BVP4C present in MATLAB
6.0.
% Starting with MATLAB 7.0, BVP4C solves multi-point BVPs
directly.
%
% The quantity of most interest is the emergent
osmolarity Os =
% 1/v(lambda). The parameters are related to another
parameter
% kappa by eta = lambda^2/(n*kappa^2). Lin and Segel
develop an
% approximate solution for Os valid for "small" n. Here
the BVP is
% solved for a range of kappa when lambda = 2 and n =
0.005. The
% computed Os is compared to the approximation of Lin and
Segel.

% Copyright 2004, The MathWorks, Inc.

Page 26
 Bvp4c_Exmamples.txt

% Known parameters, visible in nested functions.

n = 5e-2;
lambda = 2;

% Initial mesh - duplicate the interface point x = 1.

xinit = [0, 0.25, 0.5, 0.75, 1, 1, 1.25, 1.5, 1.75, 2];
lambda = 2;

sol = bvpinit(xinit,[1 1]);

fprintf(' kappa computed Os approximate Os \n')
for kappa = 2:5
eta = lambda^2/(n*kappa^2);
% After creating function handles, the new value of eta
% will be used in nested functions.
sol = bvp4c(@ex9mode,@ex9mbc,sol);

K2 = lambda*sinh(kappa/lambda)/(kappa*cosh(kappa));
approx = 1/(1 - K2);
computed = 1/sol.y(1,end);
fprintf(' %2i %10.3f %10.3f
\n',kappa,computed,approx);
end

figure
plot(sol.x,sol.y(1,:),sol.x,sol.y(2,:),'--')
legend('v(x)','C(x)')
title('A three-point BVP.')
xlabel(['\lambda = ',num2str(lambda),', \kappa =
',num2str(kappa),'.'])
ylabel('v and C')

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

function dydx = ex9mode(x,y,region)

%EX9MODE ODE function for Example 9 of the BVP tutorial.

% Here the problem is solved directly, as a three-point

Page 27
 Bvp4c_Exmamples.txt
BVP.
dydx = zeros(2,1);
dydx(1) = (y(2) - 1)/n;
% The definition of C'(x) depends on the region.
switch region
case 1 % x in [0 1]
dydx(2) = (y(1)*y(2) - x)/eta;
case 2 % x in [1 lambda]
dydx(2) = (y(1)*y(2) - 1)/eta;
end
end % ex9mode

%
-----------------------------------------------------------
---------------
function res = ex9mbc(YL,YR)
%EX9MBC Boundary conditions for Example 9 of the BVP
tutorial.
% Here the problem is solved directly, as a three-point
BVP.
res = [ YL(1,1) % v(0) = 0
YR(1,1) - YL(1,2) % continuity of v(x) at x =
1
YR(2,1) - YL(2,2) % continuity of C(x) at x =
1
YR(2,end) - 1 ]; % C(lambda) = 1
end % ex9mbc

%
-----------------------------------------------------------
---------------

end % ex9mbvp

function gasbvp
%GASBVP Exercise for Example 5 of the BVP tutorial.
% Example 8.4 of P.B. Bailey, L.F. Shampine, and P.E.
Waltman,
% Nonlinear Two Point Boundary Value Problems, Academic,
New York,
% 1968. This is a problem on an infinite interval with a
parameter
% alpha here taken to be 0.8. w(z) should decrease
Page 28
 Bvp4c_Exmamples.txt
monotonely
% from 1 to 0.
% Copyright 1999, The MathWorks, Inc.

infinity = 3;

options = bvpset('stats','on');

solinit = bvpinit(linspace(0,infinity,5),[0.5 -0.5]);

sol = bvp4c(@gasode,@gasbc,solinit,options);
z = sol.x;
w = sol.y;

clf reset
plot(z,w(1,:));
axis([0 infinity 0 1]);
xlabel('z');
ylabel('w');
title('Unsteady gas flow in a semi-infinite porous
medium.');
shg
%
-----------------------------------------------------------
---------------
function mmbvp
%MMBVP Exercise for Example 7 of the BVP tutorial.
% This is the Michaelis-Menten kinetics problem solved in
section 6.2 of
% H.B. Keller, Numerical Methods for Two-Point
Boundary-Value Problems,
% Dover, New York, 1992, for parameters epsilon = 0.1 and
k = 0.1.

% Copyright 2004, The MathWorks, Inc.

% Known parameters, visible in the nested functions.
epsilon = 0.1;
k = 0.1;
d = 0.001;
options = bvpset('stats','on');

Page 29
 Bvp4c_Exmamples.txt
solinit = bvpinit(linspace(d,1,5),[0.01 0],0.01);
sol = bvp4c(@mmode,@mmbc,solinit,options);
p = sol.parameters; % unknown parameter
xint = linspace(d,1);
Sxint = deval(sol,xint);

% Augment the solution array with the values y(0) = p,

y'(0) = 0
% to get a solution on [0, 1]. For this problem the
solution is
% flat near x = 0, but if it had been necessary for a
smooth graph,
% other values in [0,d] could have been obtained from the
series.
x = [0 xint];
y = [[p; 0] Sxint];
figure
plot(x,y(1,:))
title('Michaelis-Menten kinetics problem with coordinate
singularity.')

%
-----------------------------------------------------------
---------------
% Nested functions
%
function dydx = mmode(x,y,p)
%MMODE ODE function for the exercise of Example 7 of the
BVP tutorial.
dydx = [ y(2)
-2*(y(2)/x) + y(1)/(epsilon*(y(1) + k)) ];
end % mmode

%
-----------------------------------------------------------
---------------

function res = mmbc(ya,yb,p)

%MMBC Boundary conditions for the exercise of Example 7
of the BVP tutorial.
% The boundary conditions at x = d are that y and y'
have values yatd and
Page 30
 Bvp4c_Exmamples.txt
% ypatd obtained from series expansions. The unknown
parameter p = y(0)
% is used in the expansions. y''(d) also appears in the
expansions.
% It is evaluated as a limit in the differential
equation.
yp2atd = p /(3*epsilon*(p + k));
yatd = p + 0.5*yp2atd*d^2;
ypatd = yp2atd*d;
res = [ yb(1) - 1
ya(1) - yatd
ya(2) - ypatd ];
end % mmbc

%
-----------------------------------------------------------
---------------

end % mmbvp
function trbvp
%TRBVP Exercise for Example 3 of the BVP tutorial.
% This problem is studied in section 5.4 of B.A.
Finlayson, The Method of
% Weighted Residuals and Variational Principles,
Academic, New York, 1972.
% It arises when modelling a tubular reactor with axial
dispersion. An
% isothermal situation with n-th order irreversible
reaction leads to
%
% y'' = Pe*(y' - R*y^n)
%
% Here Pe is the axial Peclet number and R is the
reaction rate group.
% The boundary conditions are
%
% y'(0) = Pe*(y(0) - 1), y'(1) = 0.
%
% Finlayson compares results he obtains with an
orthogonal collocation
% method to results obtained by others with finite
differences when
% Pe = 1, R = 2, and n = 2. These results are compared
Page 31
 Bvp4c_Exmamples.txt
here to results
% obtained with BVP4C. BVPVAL is used to get a smoother
graph of y(x).

% Copyright 1999, The MathWorks, Inc.

% Known parameter
Pe = 1;

options = bvpset('stats','on');
solinit = bvpinit(linspace(0,1,5),[0.5 0]);
sol = bvp4c(@trode,@trbc,solinit,options,Pe);

fprintf('\n');
fprintf('Other authors report y(0) = 0.63678, y(1) =
0.45759.\n');
fprintf('Values computed are y(0) = %7.5f, y(1) =
%7.5f\n',sol.y(1,1),sol.y(1,end));

clf reset
xint = linspace(0,1);
Sxint = bvpval(sol,xint);
plot(xint,Sxint(1,:))
title('Mass transfer in a tubular reactor.')
xlabel('x')
ylabel('y')
shg
%
-----------------------------------------------------------
---------------

function dydx = trode(x,y,Pe)

%TRODE ODE function for the exercise of Example 3 of the
BVP tutorial.
dydx = [ y(2)
Pe*(y(2) + 2*y(1)^2)];

%
-----------------------------------------------------------
---------------
function res = trbc(ya,yb,Pe)
%TRBC Boundary conditions for the exercise of Example 3 of
Page 32
 Bvp4c_Exmamples.txt
the BVP tutorial.
res = [ ya(2) - Pe*(ya(1) - 1)
yb(2) ];

Page 33