Solver  Description 

Solve boundary value problems for ordinary differential equations.  
Solve boundary value problems for ordinary differential equations. 
Function  Description 

Form the initial guess for  
Evaluate the numerical solution using the output of 
An options structure contains named properties whose values
are passed to bvp4c
, and which affect problem solution.
Use these functions to create, alter, or access an options structure.
The BVP solver is designed to handle systems of ordinary differential equations
y′ = f(x, y)
where x is the independent variable, y is the dependent variable, and y′ represents the derivative of y with respect to x dy/dx.
See Choose an ODE Solver for general information about ODEs.
In a boundary value problem, the solution
of interest satisfies certain boundary conditions. These conditions
specify a relationship between the values of the solution at more
than one x. In its basic syntax, bvp4c
is designed to solve twopoint
BVPs, i.e., problems where the solution sought on an interval [a, b]
must satisfy the boundary conditions
g(y(a), y(b)) = 0
Unlike initial value problems, a boundary value problem may not have a solution, may have a finite number of solutions, or may have infinitely many solutions. As an integral part of the process of solving a BVP, you need to provide a guess for the required solution. The quality of this guess can be critical for the solver performance and even for a successful computation.
There may be other difficulties when solving BVPs, such as problems imposed on infinite intervals or problems that involve singular coefficients. Often BVPs involve unknown parameters p that have to be determined as part of solving the problem
y′ = f(x, y, p)
g(y(a), y(b), p) = 0
In this case, the boundary conditions must suffice to determine the value of p.
The function bvp4c
solves
twopoint boundary value problems for ordinary differential equations
(ODEs). It integrates a system of firstorder ordinary differential
equations
y′ = f(x, y)
on the interval [a, b], subject to general twopoint boundary conditions
bc(y(a),y(b)) = 0
It can also accommodate other types of BVP problems, such as those that have any of the following:
Unknown parameters
Singularities in the solutions
Multipoint conditions
In this case, the number of boundary conditions must be sufficient to determine the solution and the unknown parameters.
bvp4c
produces a solution that is continuous
on [a,b] and has a continuous
first derivative there. You can use the function deval
and the output of bvp4c
to
evaluate the solution at specific points on the interval of integration.
bvp4c
is a finite difference code that implements
the 3stage Lobatto IIIa formula. This is a collocation formula and
the collocation polynomial provides a C^{1}continuous
solution that is fourthorder accurate uniformly in the interval of
integration. Mesh selection and error control are based on the residual
of the continuous solution.
The collocation technique uses a mesh of points to divide the interval of integration into subintervals. The solver determines a numerical solution by solving a global system of algebraic equations resulting from the boundary conditions, and the collocation conditions imposed on all the subintervals. The solver then estimates the error of the numerical solution on each subinterval. If the solution does not satisfy the tolerance criteria, the solver adapts the mesh and repeats the process. The user must provide the points of the initial mesh as well as an initial approximation of the solution at the mesh points.
The basic syntax of the BVP solver is
sol = bvp4c(odefun,bcfun,solinit)
The input arguments are
 A function handle that evaluates the differential equations. It has the basic form dydx = odefun(x,y) where  
 Handle to a function that evaluates the residual in the boundary conditions. It has the basic form res = bcfun(ya,yb) where  
 Structure with fields  
 Ordered nodes of the initial mesh. Boundary conditions
are imposed at a =  
 Initial guess for the solution with  
The structure can have any
name, but the fields must be named 
The output argument sol
is a structure created
by the solver. In the basic case the structure has fields x
, y
, yp
,
and solver
.
 Nodes of the mesh selected by 
 Approximation to y(x)
at the mesh points of 
 Approximation to y′(x)
at the mesh points of 


The structure sol
returned by bvp4c
contains
an additional field if the problem involves unknown parameters:
 Value of unknown parameters, if present, found by the solver. 
The function deval
uses
the output structure sol
to evaluate the numerical
solution at any point from [a,b]
.
For more advanced applications, you can specify solver options
by passing an input argument options
.
The default integration properties in the BVP solver bvp4c
are
selected to handle common problems. In some cases, you can improve
solver performance by overriding these defaults. You do this by supplying bvp4c
with
an options
structure that specifies one or more
property values.
For example, to change the value of the relative error tolerance
of bvp4c
from the default value of 1e3
to 1e4
,
Create an options structure using the function bvpset
by entering
options = bvpset('RelTol', 1e4);
Pass the options
structure
to bvp4c
as follows:
sol = bvp4c(odefun,bcfun,solinit,options)
For a complete description of the available options, see the
reference page for bvpset
.
Solving the Problem. This example determines the fourth eigenvalue of Mathieu's Equation.
It illustrates how to write secondorder differential equations as
a system of two firstorder ODEs and how to use bvp4c
to
determine an unknown parameter λ.
The task is to compute the fourth (q = 5) eigenvalue lambda λ of Mathieu's equation
y′′ + (λ – 2 q cos 2x)y = 0
Because the unknown parameter λ is present, this secondorder differential equation is subject to three boundary conditions
y(0) = 1
y′(0) = 0
y′(π) = 0
Note:
The file, 
Rewrite the problem
as a firstorder system. To use bvp4c
,
you must rewrite the equations as an equivalent system of firstorder
differential equations. Using a substitution y_{1} = y and y_{2} = y′,
the differential equation is written as a system of two firstorder
equations
y_{1}′ = y_{2}
y_{2}′ = –(λ – 2q cos 2x)y_{1}
Note that the differential equations depend on the unknown parameter λ. The boundary conditions become
y_{1}(0) – 1 = 0
y_{2}(0) = 0
y_{2}(π) = 0
Code the system of
firstorder ODEs. Once you represent the equation as a
firstorder system, you can code it as a function that bvp4c
can
use. Because there is an unknown parameter, the function must be of
the form
dydx = odefun(x,y,parameters)
The following code represents the system in the function, mat4ode
.
Variable q
is shared with the outer function:
function dydx = mat4ode(x,y,lambda) dydx = [ y(2) (lambda  2*q*cos(2*x))*y(1) ]; end % End nested function mat4ode
Code the boundary conditions function. You must also code the boundary conditions in a function. Because there is an unknown parameter, the function must be of the form
res = bcfun(ya,yb,parameters)
The code below represents the boundary conditions in the function, mat4bc
.
function res = mat4bc(ya,yb,lambda) res = [ ya(2) yb(2) ya(1)1 ];
Create an initial guess. To form the guess structure solinit
with bvpinit
, you need to provide initial
guesses for both the solution and the unknown parameter.
The function mat4init
provides an initial
guess for the solution. mat4init
uses y =
cos4x because this function satisfies the boundary
conditions and has the correct qualitative behavior (the correct number
of sign changes).
function yinit = mat4init(x) yinit = [ cos(4*x) 4*sin(4*x) ];
In the call to bvpinit
, the third argument, lambda
,
provides an initial guess for the unknown parameter λ.
lambda = 15; solinit = bvpinit(linspace(0,pi,10),@mat4init,lambda);
This example uses the @
symbol to pass mat4init
as
a function handle to bvpinit
.
Apply the BVP solver. The mat4bvp
example
calls bvp4c
with the functions mat4ode
and mat4bc
and
the structure solinit
created with bvpinit
.
sol = bvp4c(@mat4ode,@mat4bc,solinit);
View the results. Complete the example by displaying the results:
Print the value of the unknown parameter λ found
by bvp4c
.
fprintf('Fourth eigenvalue is approximately %7.3f.\n',... sol.parameters)
Use deval
to
evaluate the numerical solution at 100 equally spaced points in the
interval [0, π], and plot its first component. This component
approximates y(x).
xint = linspace(0,pi); Sxint = deval(sol,xint); plot(xint,Sxint(1,:)) axis([0 pi 1 1.1]) title('Eigenfunction of Mathieu''s equation.') xlabel('x') ylabel('solution y')
The following plot shows the eigenfunction associated with the final eigenvalue λ = 17.097.
Finding Unknown Parameters. The bvp4c
solver can
find unknown parameters p for problems of the form
y′ = f(x,y,p)
bc(y(a), y (b),p) = 0
You must provide bvp4c
an initial guess for
any unknown parameters in the vector solinit.parameters
.
When you call bvpinit
to create
the structure solinit
, specify the initial guess
as a vector in the additional argument parameters
.
solinit = bvpinit(x,v,parameters)
The bvp4c
function arguments odefun
and bcfun
must
each have a third argument.
dydx = odefun(x,y,parameters) res = bcfun(ya,yb,parameters)
While solving the differential equations, bvp4c
adjusts
the value of unknown parameters to satisfy the boundary conditions.
The solver returns the final values of these unknown parameters in sol.parameters
.
Evaluating the Solution. The collocation method implemented in bvp4c
produces
a C^{1}continuous solution over the whole
interval of integration [a,b].
You can evaluate the approximate solution, S(x),
at any point in [a,b] using
the helper function deval
and the structure sol
returned
by bvp4c
.
Sxint = deval(sol,xint)
The deval
function is vectorized. For a
vector xint
, the i
th column
of Sxint
approximates the solution y(xint(i)).
Introduction. To solve a boundary value problem, you need to provide an initial guess for the solution. The quality of your initial guess can be critical to the solver performance, and to being able to solve the problem at all. However, coming up with a sufficiently good guess can be the most challenging part of solving a boundary value problem. Certainly, you should apply the knowledge of the problem's physical origin. Often a problem can be solved as a sequence of relatively simpler problems, i.e., a continuation.
This example shows how to use continuation to:
Solve a difficult BVP
Verify a solution's consistent behavior
Using Continuation to Solve a BVP. This example solves the differential equation
εy′′ + xy′ = επ^{2} cos(πx) – πx sin(πx)
for ε = 10^{–4}, on the interval [–1 1], with boundary conditions y(–1) = –2 and y(1) = 0. For 0 < ε <1, the solution has a transition layer at x = 0. Because of this rapid change in the solution for small values of ε, the problem becomes difficult to solve numerically.
The example solves the problem as a sequence of relatively simpler problems, i.e., a continuation. The solution of one problem is used as the initial guess for solving the next problem.
Note:
The file, 
Note: This problem appears in [1] to illustrate the mesh selection capability of a well established BVP code COLSYS. 
Code the ODE and
boundary condition functions. Code the differential equation
and the boundary conditions as functions that bvp4c
can
use:
The code below represents the differential equation and the
boundary conditions in the functions shockODE
and shockBC
.
Note that shockODE
is vectorized to improve solver
performance. The additional parameter ε is represented by e
and
is shared with the outer function.
function dydx = shockODE(x,y) pix = pi*x; dydx = [ y(2,:) x/e.*y(2,:)  pi^2*cos(pix)  pix/e.*sin(pix) ]; end % End nested function shockODE function res = shockBC(ya,yb) res = [ ya(1)+2 yb(1) ]; end % End nested function shockBC
Provide analytical
partial derivatives. For this problem, the solver benefits
from using analytical partial derivatives. The code below represents
the derivatives in functions shockJac
and shockBCJac
.
function jac = shockJac(x,y) jac = [ 0 1 0 x/e ]; end % End nested function shockJac function [dBCdya,dBCdyb] = shockBCJac(ya,yb) dBCdya = [ 1 0 0 0 ]; dBCdyb = [ 0 0 1 0 ]; end % End nested function shockBCJac
shockJac shares e
with the outer function.
Tell bvp4c
to use these functions to evaluate
the partial derivatives by setting the options FJacobian and BCJacobian.
Also set 'Vectorized'
to 'on'
to
indicate that the differential equation function shockODE
is
vectorized.
options = bvpset('FJacobian',@shockJac,... 'BCJacobian',@shockBCJac,... 'Vectorized','on');
Create an initial
guess. You must provide bvp4c
with a
guess structure that contains an initial mesh and a guess for values
of the solution at the mesh points. A constant guess of y(x)
≡ 1 and y′(x)
≡ 0, and a mesh of five equally spaced points on [–1
1] suffice to solve the problem for ε = 10^{–2}.
Use bvpinit
to form the guess
structure.
sol = bvpinit([1 0.5 0 0.5 1],[1 0]);
Use continuation
to solve the problem. To obtain the solution for the parameter
ε = 10^{–4}, the example uses
continuation by solving a sequence of problems for ε = 10^{–2},
10^{–3}, 10^{–4}.
The solver bvp4c
does not perform continuation
automatically, but the code's user interface has been designed to
make continuation easy. The code uses the output sol
that bvp4c
produces
for one value of e
as the guess in the next iteration.
e = 0.1; for i=2:4 e = e/10; sol = bvp4c(@shockODE,@shockBC,sol,options); end
View the results. Complete the example by displaying the final solution
plot(sol.x,sol.y(1,:)) axis([1 1 2.2 2.2]) title(['There is a shock at x = 0 when \epsilon = '... sprintf('%.e',e) '.']) xlabel('x') ylabel('solution y')
Using Continuation to Verify Consistency. FalknerSkan BVPs arise from similarity solutions of viscous, incompressible, laminar flow over a flat plate. An example is
f ′′′ + ff ′′ + β(1 – (f ′)^{2}) = 0
for β = 0.5 on the interval [0, ∞] with boundary conditions f(0) = 0, f ′(0) = 0, and f ′(∞) = 1.
The BVP cannot be solved on an infinite interval, and it would be impractical to solve it for even a very large finite interval. So, the example tries to solve a sequence of problems posed on increasingly larger intervals to verify the solution's consistent behavior as the boundary approaches ∞.
The example imposes the infinite boundary condition at a finite
point called infinity
. The example then uses continuation
in this end point to get convergence for increasingly larger values
of infinity
. It uses bvpinit
to
extrapolate the solution sol
for one value of infinity
as
an initial guess for the new value of infinity
.
The plot of each successive solution is superimposed over those of
previous solutions so they can easily be compared for consistency.
Note:
The file, 
Code the ODE and
boundary condition functions. Code the differential equation
and the boundary conditions as functions that bvp4c
can
use. The problem parameter beta
is shared with
the outer function.
function dfdeta = fsode(eta,f) dfdeta = [ f(2) f(3) f(1)*f(3)  beta*(1  f(2)^2) ]; end % End nested function fsode function res = fsbc(f0,finf) res = [f0(1) f0(2) finf(2)  1]; end % End nested function fsbc
Create an initial
guess. You must provide bvp4c
with a
guess structure that contains an initial mesh and a guess for values
of the solution at the mesh points. A crude mesh of five points and
a constant guess that satisfies the boundary conditions are good enough
to get convergence when infinity = 3
.
infinity = 3; maxinfinity = 6; solinit = bvpinit(linspace(0,infinity,5),[0 0 1]);
Solve on the initial
interval. The example obtains the solution for infinity
= 3
. It then prints the computed value of f ′′(0) for comparison with the value
reported by Cebeci and Keller [2]:
sol = bvp4c(@fsode,@fsbc,solinit); eta = sol.x; f = sol.y; fprintf('\n'); fprintf('Cebeci & Keller report that f''''(0) = 0.92768.\n') fprintf('Value computed using infinity = %g is %7.5f.\n', ... infinity,f(3,1))
The example prints
Cebeci & Keller report that f''(0) = 0.92768. Value computed using infinity = 3 is 0.92915.
Setup the figure and plot the initial solution.
figure plot(eta,f(2,:),eta(end),f(2,end),'o'); axis([0 maxinfinity 0 1.4]); title('FalknerSkan equation, positive wall shear, \beta = 0.5.') xlabel('\eta') ylabel('df/d\eta') hold on drawnow shg
Use continuation
to solve the problem and plot subsequent solutions. The
example then solves the problem for infinity = 4
, 5
, 6
.
It uses bvpinit
to extrapolate
the solution sol
for one value of infinity
as
an initial guess for the next value of infinity
.
For each iteration, the example prints the computed value of f ′′(0) and superimposes a plot of
the solution in the existing figure.
for Bnew = infinity+1:maxinfinity solinit = bvpinit(sol,[0 Bnew]); % Extend solution to Bnew. sol = bvp4c(@fsode,@fsbc,solinit); eta = sol.x; f = sol.y; fprintf('Value computed using infinity = %g is %7.5f.\n', ... Bnew,f(3,1)) plot(eta,f(2,:),eta(end),f(2,end),'o'); drawnow end hold off
The example prints
Value computed using infinity = 4 is 0.92774. Value computed using infinity = 5 is 0.92770. Value computed using infinity = 6 is 0.92770.
Note that the values approach 0.92768 as reported by Cebeci and Keller. The superimposed plots confirm the consistency of the solution's behavior.
Introduction. The function bvp4c
solves
a class of singular BVPs of the form
$$\begin{array}{c}{y}^{\prime}=\frac{1}{x}Sy+f(x,y)\\ 0=g\left(y(0),y(b)\right)\end{array}$$  (121) 
It can also accommodate unknown parameters for problems of the form
$$\begin{array}{c}{y}^{\prime}=\frac{1}{x}Sy+f(x,y,p)\\ 0=\left(g(y(0),y(b),p\right)\end{array}$$
Singular problems must be posed on an interval [0,b]
with b > 0. Use bvpset
to
pass the constant matrix S to bvp4c
as
the value of the 'SingularTerm'
integration property.
Boundary conditions at x = 0 must be consistent
with the necessary condition for a smooth solution, Sy(0)
= 0. An initial guess should also satisfy this necessary condition.
When you solve a singular BVP using
sol = bvp4c(@odefun,@bcfun,solinit,options)
bvp4c
requires that your function odefun(x,y)
return
only the value of the f(x, y)
term in Equation 52.
Emden's equation. Emden's equation arises in modeling a spherical body of gas. The PDE of the model is reduced by symmetry to the ODE
$${{y}^{\prime}}^{\prime}+\frac{2}{x}{y}^{\prime}+{y}^{5}=0$$
on an interval [0,1]. The coefficient 2/x is singular at x = 0, but symmetry implies the boundary condition y′(0) = 0. With this boundary condition, the term
$$\frac{2}{x}{y}^{\prime}(0)$$
is welldefined as x approaches 0. For the boundary condition $$y(1)=\sqrt{3/2}$$, this BVP has the analytical solution
$$y(x)={\left(1+\frac{{x}^{2}}{3}\right)}^{1/2}$$
Note:
The file, 
Rewrite the problem as a firstorder system and identify the singular term. Using a substitution y_{1} = y and y_{2} = y′, write the differential equation as a system of two firstorder equations
$$\begin{array}{l}{y}_{1}{}^{\prime}={y}_{2}\\ {y}_{2}{}^{\prime}=\frac{2}{x}{y}_{2}{y}_{1}^{5}\end{array}$$
The boundary conditions become
$$\begin{array}{l}{y}_{2}(0)=0\\ {y}_{1}(1)=\sqrt{3/2}\end{array}$$
Writing the ODE system in a vectormatrix form
$$\left[\begin{array}{c}{y}_{1}{}^{\prime}\\ {y}_{2}{}^{\prime}\end{array}\right]=\frac{1}{x}\left[\begin{array}{cc}0& 0\\ 0& 2\end{array}\right]\left[\begin{array}{c}{y}_{1}\\ {y}_{2}\end{array}\right]+\left[\begin{array}{c}{y}_{2}\\ {y}_{1}^{5}\end{array}\right]$$
the terms of Equation 52 are identified as
$$S=\left[\begin{array}{cc}0& 0\\ 0& 2\end{array}\right]$$
and
$$f(x,y)=\left[\begin{array}{c}{y}_{2}\\ {y}_{1}^{5}\end{array}\right]$$
Code the ODE and
boundary condition functions. Code the differential equation
and the boundary conditions as functions that bvp4c
can
use.
function dydx = emdenode(x,y) dydx = [ y(2) y(1)^5 ]; function res = emdenbc(ya,yb) res = [ ya(2) yb(1)  sqrt(3)/2 ];
Setup integration
properties. Use the matrix as the value of the 'SingularTerm'
integration
property.
S = [0,0;0,2]; options = bvpset('SingularTerm',S);
Create an initial guess. This example starts with a mesh of five points and a constant guess for the solution.
$$\begin{array}{l}{y}_{1}(x)\equiv \sqrt{3/2}\\ {y}_{2}(x)\equiv 0\end{array}$$
Use bvpinit
to form the
guess structure
guess = [sqrt(3)/2;0]; solinit = bvpinit(linspace(0,1,5),guess);
Solve the problem. Use
the standard bvp4c
syntax to solve the problem.
sol = bvp4c(@emdenode,@emdenbc,solinit,options);
View the results. This problem has an analytical solution
$$y(x)={\left(1+\frac{{x}^{2}}{3}\right)}^{1/2}$$
The example evaluates the analytical solution at 100 equally
spaced points and plots it along with the numerical solution computed
using bvp4c
.
x = linspace(0,1); truy = 1 ./ sqrt(1 + (x.^2)/3); plot(x,truy,sol.x,sol.y(1,:),'ro'); title('Emden problem  BVP with singular term.') legend('Analytical','Computed'); xlabel('x'); ylabel('solution y');
In multipoint boundary value problems, the solution of interest
satisfies conditions at points inside the interval of integration.
The bvp4c
function is useful
in solving such problems.
The following example shows how the multipoint capability in bvp4c
can
improve efficiency when you are solving a nonsmooth problem. The following
equations are solved on 0 ≤ x ≤ λ for constant
parameters n
, κ, λ > 1, and η
= λ^{2}/(n × κ^{2}).
These are subject to boundary conditions v(0) = 0
and
C(λ) = 1:
v' = (C  1)/n
C' = (v * C  min(x,1))/$$\eta $$
The term min(x,1)
is not smooth at xc
= 1
, and this can affect the solver's efficiency. By introducing
an interface point at xc = 1
, smooth solutions
can be obtained on [0,1] and [1,λ]. To get a continuous solution
over the entire interval [0,λ], the example imposes matching
conditions at the interface.
Note:
The file, 
The example takes you through the following steps:
Determine the interfaces and divide
the interval of integration into regions. Introducing an
interface point at xc = 1
divides the problem into
two regions in which the solutions remain smooth. The differential
equations for the two regions are
Region 1: 0 ≤ x ≤ 1
v' = (C  1)/n
C' = (v * C  x)/$$\eta $$
Region 2: 1 ≤ x ≤ λ
v' = (C  1)/n
C' = (v * C  1)/$$\eta $$
Note that the interface xc = 1
is included
in both regions. At xc = 1
, bvp4c
produces
a left and right solution.
These solutions are denoted as v(1)
, C(1)
and v(1+)
, C(1+)
respectively.
Determine the boundary conditions. Solving
two firstorder differential equations in two regions requires imposing
four boundary conditions. Two of these conditions come from the original
formulation; the others enforce the continuity of the solution across
the interface xc = 1
:
v(0) = 0
C($$\lambda $$)  1 = 0
v(1)  v(1+) = 0
C(1)  C(1+) = 0
Here, v(1), C(1)
and v(1+), C(1+)
denote
the left and right solution at the interface.
Code the derivative function. In
the derivative function, y(1)
corresponds to v(x)
,
and y(2)
corresponds to C(x)
.
The additional input argument region
identifies
the region in which the derivative is evaluated. bvp4c
enumerates
regions from left to right, starting with 1. Note that the problem
parameters n
and η are shared with the outer
function:
function dydx = f(x,y,region) 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 % End nested function f
Code the boundary conditions function. For
multipoint BVPs, the arguments of the boundary conditions function, YL
and YR
,
become matrices. In particular, the k
th column YL(:,k)
represents
the solution at the left boundary of the k
th region.
Similarly, YR(:,k)
represents the solution at the
right boundary of the k
th region.
In the example, y(0) is approximated by YL(:,1)
,
while y(λ) is approximated by YR(:,end)
.
Continuity of the solution at the internal interface requires that YR(:,1)
= YL(:,2)
. Nested function bc
computes
the residual in the boundary conditions:
function res = bc(YL,YR)
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 % End nested function bc
Create an initial guess. For
multipoint BVPs, when creating an initial guess using bvpinit
,
use double entries in xinit
for the interface point xc
.
This example uses a constant guess yinit = [1;1]
:
xc = 1; xinit = [0, 0.25, 0.5, 0.75, xc, xc, 1.25, 1.5, 1.75, 2]; solinit = bvpinit(xinit,yinit)
For multipoint BVPs, you can use different guesses in different
regions. To do that, you specify the initial guess for y
as
a function using the following syntax:
solinit = bvpinit(xinit,@yinitfcn)
The initial guess function must have the following general form:
function y = yinitfcn(x,region) switch region case 1 % x in [0, 1] y = [1;1]; % initial guess for y(x) $$0\le \text{x}\le \text{1}$$ case 2 % x in [1, $$\lambda $$] y = [1;1]; % initial guess for y(x), $$1\le \text{x}\le \lambda $$ end
Apply the solver. The bvp4c
function uses the same syntax for
multipoint BVPs as it does for twopoint BVPs:
sol = bvp4c(@f,@bc,solinit);
The mesh points returned in sol.x
are adapted
to the solution behavior, but the mesh still includes a double entry
for the interface point xc = 1
. Corresponding columns
of sol.y
represent the left and right solution
at xc
.
View the results. Using deval
, the solution can be evaluated
at any point in the interval of integration.
Note that, with the left and right values computed at the interface,
the solution is not uniquely defined at xc = 1
.
When evaluating the solution exactly at the interface, deval
issues
a warning and returns the average of the left and right solution values.
Call deval
at xceps(xc)
and xc+eps(xc)
to
get the limit values at xc
.
The example plots the solution approximated at the mesh points selected by the solver:
plot(sol.x,sol.y(1,:),sol.x,sol.y(2,:),'') legend('v(x)','C(x)') title(['A threepoint BVP solved with BVP4C']) xlabel(['\lambda = 2, \kappa = 5.']) ylabel('v and C')
The following additional examples are available. Type
edit examplename
to view the code and
examplename
to run the example.
Example Name  Description 

 Emden's equation, a singular BVP 
 FalknerSkan BVP on an infinite interval 
 Fourth eigenfunction of Mathieu's equation 
 Solution with a shock layer near x = 0 
 BVP with exactly two solutions 
 Threepoint boundary value problem 
Additional examples are provided by "Tutorial on Solving
BVPs with BVP4C," available at http://www.mathworks.com/bvp_tutorial
.
[1] Ascher, U., R. Mattheij, and R. Russell, "Numerical Solution of Boundary Value Problems for Ordinary Differential Equations," SIAM, Philadelphia, PA, 1995, p. 372.
[2] Cebeci, T. and H. B. Keller, "Shooting and Parallel Shooting Methods for Solving the FalknerSkan Boundarylayer Equation," J. Comp. Phys., Vol. 7, 1971, pp. 289300.