Solving Boundary Value Problems
In a boundary value problem (BVP), the goal is to find a solution to an ordinary differential equation (ODE) that also satisfies certain specified boundary conditions. The boundary conditions specify a relationship between the values of the solution at two or more locations in the interval of integration. In the simplest case, the boundary conditions apply at the beginning and end (or boundaries) of the interval.
The MATLAB® BVP solvers
designed to handle systems of ODEs of the form
x is the independent variable
y is the dependent variable
y′ represents the derivative of y with respect to x, also written as dy/dx
In the simplest case of a two-point BVP, the solution to the ODE is sought on an interval [a, b], and must satisfy the boundary conditions
To specify the boundary conditions for a given BVP, you must:
Write a function of the form
res = bcfun(ya,yb), or use the form
res = bcfun(ya,yb,p)if there are unknown parameters involved. You supply this function to the solver as the second input argument. The function returns
res, which is the residual value of the solution at the boundary point. For example, if y(a) = 1 and y(b) = 0, then the boundary condition function is
function res = bcfun(ya,yb) res = [ya(1)-1 yb(1)]; end
In the initial guess for the solution, the first and last points in the mesh specify the points at which the boundary conditions are enforced. For the above boundary conditions, you can specify
bvpinit(linspace(a,b,5),yinit)to enforce the boundary conditions at a and b.
The BVP solvers in MATLAB also can accommodate other types of problems that have:
Unknown parameters p
Singularities in the solutions
Multipoint conditions (internal boundaries that separate the integration interval into several regions)
In the case of multipoint boundary conditions, the boundary
conditions apply at more than two points in the interval of integration. For
example, the solution might be required to be zero at the beginning, middle, and end
of the interval. See
bvpinit for details on how to
specify multiple boundary conditions.
Initial Guess of Solution
Unlike initial value problems, a boundary value problem can have:
A finite number of solutions
Infinitely many solutions
An important part of the process of solving a BVP is providing a guess for the required solution. The quality of this guess can be critical for the solver performance and even for a successful computation.
bvpinit function to create a
structure for the initial guess of the solution. The solvers
bvp5c accept this structure
as the third input argument.
Creating a good initial guess for the solution is more an art than a science. However, some general guidelines include:
Have the initial guess satisfy the boundary conditions, since the solution is required to satisfy them as well. If the problem contains unknown parameters, then the initial guess for the parameters also should satisfy the boundary conditions.
Try to incorporate as much information about the physical problem or expected solution into the initial guess as possible. For example, if the solution is supposed to oscillate or have a certain number of sign changes, then the initial guess should as well.
Consider the placement of the mesh points (the x-coordinates of the initial guess of the solution). The BVP solvers adapt these points during the solution process, so you do not need to specify too many mesh points. Best practice is to specify a few mesh points placed near where the solution changes rapidly.
If there is a known, simpler solution on a smaller interval, then use it as an initial guess on a larger interval. Often you can solve a problem as a series of relatively simpler problems, a practice called continuation. With continuation, a series of simple problems are connected by using the solution of one problem as the initial guess to solve the next problem.
Finding Unknown Parameters
Often BVPs involve unknown parameters p that have to be determined as part of solving the problem. The ODE and boundary conditions become
In this case, the boundary conditions must suffice to determine the values of the parameters p.
You must provide the solver with an initial guess for any unknown parameters. When
bvpinit to create the structure
solinit, specify the initial guess as a vector in the third
solinit = bvpinit(x,v,parameters)
Additionally, the functions
that encode the ODE equations and boundary conditions must each have a third
dydx = odefun(x,y,parameters) res = bcfun(ya,yb,parameters)
While solving the differential equations, the solver adjusts the value of the
unknown parameters to satisfy the boundary conditions. The solver returns the final
values of these unknown parameters in
bvp5c can solve a class of
singular BVPs of the form
The solvers can also accommodate unknown parameters for problems of the form
Singular problems must be posed on an interval [0,b] with
b > 0. Use
bvpset to pass the constant matrix
S to the solver as the value of the
'SingularTerm' option. Boundary conditions at
x = 0 must be consistent with the necessary condition for a
smooth solution, Sy(0) = 0. The initial guess of the solution
also should satisfy this condition.
When you solve a singular BVP, the solver requires that your function
odefun(x,y) return only the value of the
f(x, y) term in the
equation. The term involving S is handled by the solver
separately using the
BVP Solver Selection
MATLAB includes the solvers
bvp5c to solve BVPs. In most cases you can use the solvers
interchangeably. The main difference between the solvers is that
bvp4c implements a fourth-order formula, while
bvp5c implements a fifth-order formula.
bvp5c function is used exactly
bvp4c, with the exception of the meaning of error
tolerances between the two solvers. If S(x)
approximates the solution y(x),
controls the residual
|S′(x) – f(x,S(x))|.
This approach indirectly controls the true error
control the true error directly.
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, then the solver adapts the mesh and repeats the process. You must provide the points of the initial mesh, as well as an initial approximation of the solution at the mesh points.
bvp5c is a finite difference code that
implements the four-stage Lobatto IIIa formula. This is a
collocation formula and the collocation polynomial provides a
C1-continuous solution that is
fifth-order accurate uniformly in [a,b]. The formula is implemented
as an implicit Runge-Kutta formula.
solves the algebraic equations directly, whereas
bvp4c uses analytical condensation.
bvp4c handles unknown parameters directly,
bvp5c augments the system with trivial
differential equations for unknown parameters.
Evaluating the Solution
The collocation methods implemented in
bvp5c produce C1-continuous
solutions over the 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
the solver. For example, to evaluate the solution
sol at the mesh
xint, use the command
Sxint = deval(sol,xint)
deval function is vectorized. For a vector
ith column of
Sxint approximates the solution
BVP Examples and Files
Several available example files serve as excellent starting points for most common BVP problems. To easily explore and run examples, simply use the Differential Equations Examples app. To run this app, type
This table contains a list of the available BVP example files, as well as the solvers and the options they use.
Emden's equation, a singular BVP
Falkner-Skan BVP on an infinite interval
Fourth eigenfunction of Mathieu's equation
Example comparing the errors controlled by
Compare bvp4c and bvp5c Solvers (
Compare bvp4c and bvp5c Solvers (
Solution with a shock layer near x = 0
BVP with exactly two solutions
Three-point boundary value problem
 Ascher, U., R. Mattheij, and R. Russell. “Numerical Solution of Boundary Value Problems for Ordinary Differential Equations.” Philadelphia, PA: SIAM, 1995, p. 372.
 Shampine, L.F., and J. Kierzenka. "A BVP Solver based on residual control and the MATLAB PSE." ACM Trans. Math. Softw. Vol. 27, Number 3, 2001, pp. 299–316.
 Shampine, L.F., M.W. Reichelt, and J. Kierzenka. "Solving Boundary Value Problems for Ordinary Differential Equations in MATLAB with bvp4c." MATLAB File Exchange, 2004.
 Shampine, L.F., and J. Kierzenka. "A BVP Solver that Controls Residual and Error." J. Numer. Anal. Ind. Appl. Math. Vol. 3(1-2), 2008, pp. 27–41.