This is the MATLAB^{®} PDE solver.
PDE InitialBoundary Value Problem Solver  Description 

Solve initialboundary value problems for systems of parabolic and elliptic PDEs in one space variable and time. 
PDE Helper Function  Description 

Evaluate the numerical solution of a PDE using the output
of 
pdepe
solves systems of parabolic and elliptic
PDEs in one spatial variable x and time t,
of the form
$$c\left(x,t,u,\frac{\partial u}{\partial x}\right)\frac{\partial u}{\partial t}={x}^{m}\frac{\partial}{\partial x}\left({x}^{m}f\left(x,t,u,\frac{\partial u}{\partial x}\right)\right)+s\left(x,t,u,\frac{\partial u}{\partial x}\right).$$  (115) 
The PDEs hold for t_{0} ≤ t ≤ t_{f} and a ≤ x ≤ b. The interval [a, b] must be finite. m can be 0, 1, or 2, corresponding to slab, cylindrical, or spherical symmetry, respectively. If m > 0, then a ≥ 0 must also hold.
In Equation 115, f(x,t,u,∂u/∂x) is a flux term and s(x,t,u,∂u/∂x) is a source term. The flux term must depend on ∂u/∂x. The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix c(x,t,u,∂u/∂x). The diagonal elements of this matrix are either identically zero or positive. An element that is identically zero corresponds to an elliptic equation and otherwise to a parabolic equation. There must be at least one parabolic equation. An element of c that corresponds to a parabolic equation can vanish at isolated values of x if they are mesh points. Discontinuities in c and/or s due to material interfaces are permitted provided that a mesh point is placed at each interface.
At the initial time t = t_{0}, for all x the solution components satisfy initial conditions of the form
$$u(x,{t}_{0})={u}_{0}(x).$$  (116) 
At the boundary x = a or x = b, for all t the solution components satisfy a boundary condition of the form
$$p(x,t,u)+q(x,t)f\left(x,t,u,\frac{\partial u}{\partial x}\right)=0.$$  (117) 
q(x,t) is a diagonal matrix with elements that are either identically zero or never zero. Note that the boundary conditions are expressed in terms of the f rather than partial derivative of u with respect to x ∂u/∂x. Also, of the two coefficients, only p can depend on u.
The MATLAB PDE solver, pdepe
,
solves initialboundary value problems for systems of parabolic and
elliptic PDEs in the one space variable x and time t.
There must be at least one parabolic equation in the system.
The pdepe
solver converts the PDEs to ODEs
using a secondorder accurate spatial discretization based on a fixed
set of userspecified nodes. The discretization method is described
in [9]. The
time integration is done with ode15s
.
The pdepe
solver exploits the capabilities of ode15s
for
solving the differentialalgebraic equations that arise when Equation 115 contains
elliptic equations, and for handling Jacobians with a specified sparsity
pattern. ode15s
changes both the time step and
the formula dynamically.
After discretization, elliptic equations give rise to algebraic
equations. If the elements of the initial conditions vector that correspond
to elliptic equations are not "consistent" with the
discretization, pdepe
tries to adjust them before
beginning the time integration. For this reason, the solution returned
for the initial time may have a discretization error comparable to
that at any other time. If the mesh is sufficiently fine, pdepe
can
find consistent initial conditions close to the given ones. If pdepe
displays
a message that it has difficulty finding consistent initial conditions,
try refining the mesh. No adjustment is necessary for elements of
the initial conditions vector that correspond to parabolic equations.
The basic syntax of the solver is:
sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan)
Note: Correspondences given are to terms used in Initial Value Problems. 
The input arguments are
 Specifies the symmetry of the problem. 
 Function that defines the components of the PDE. It computes the terms c, f, and s in Equation 115, and has the form [c,f,s] = pdefun(x,t,u,dudx) where 
 Function that evaluates the initial conditions. It has the form u = icfun(x) When called with an argument 
 Function that evaluates the terms p and q of the boundary conditions. It has the form [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t) where 
 Vector [ Secondorder
approximation to the solution is made on the mesh specified in The elements of 
 Vector [
The elements of 
The output argument sol
is a threedimensional
array, such that
sol
(:
,:
,k
)
approximates component k
of the solution u.
sol
(i
,:
,k
)
approximates component k
of the solution at time tspan
(i
)
and mesh points xmesh(:)
.
sol
(i,j,k
) approximates
component k
of the solution at time tspan
(i
)
and the mesh point xmesh(j)
.
For more advanced applications, you can also specify as input arguments solver options and additional parameters that are passed to the PDE functions.
 Structure of optional parameters that change the default integration properties. This is the seventh input argument. sol = pdepe(m,pdefun,icfun,bcfun,... xmesh,tspan,options) See Integrator Options for more information. 
The default integration properties in the MATLAB PDE solver
are selected to handle common problems. In some cases, you can improve
solver performance by overriding these defaults. You do this by supplying pdepe
with
one or more property values in an options structure.
sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)
Use odeset
to create
the options structure. Only those options of the underlying ODE solver
shown in the following table are available for pdepe
.
The defaults obtained by leaving off the input argument options
are
generally satisfactory. Integrator Options tells you how to create
the structure and describes the properties.
PDE Properties
Properties Category  Property Name 

Error control 

Stepsize 

Solving the Equation. This example illustrates the straightforward formulation, solution, and plotting of the solution of a single PDE
$${\pi}^{2}\frac{\partial u}{\partial t}=\frac{{\partial}^{2}u}{\partial {x}^{2}}.$$
This equation holds on an interval 0 ≤ x ≤ 1 for times t ≥ 0. At t = 0, the solution satisfies the initial condition
$$u(x,0)=\mathrm{sin}\pi x.$$
At x = 0 and x = 1, the solution satisfies the boundary conditions
$$\begin{array}{l}u(0,t)=0\\ \pi {e}^{t}+\frac{\partial u}{\partial x}(1,t)=0.\end{array}$$
Note:
The file, 
Rewrite the PDE. Write the PDE in the form
$$c\left(x,t,u,\frac{\partial u}{\partial x}\right)\frac{\partial u}{\partial t}={x}^{m}\frac{\partial}{\partial x}\left({x}^{m}f\left(x,t,u,\frac{\partial u}{\partial x}\right)\right)+s\left(x,t,u,\frac{\partial u}{\partial x}\right).$$
This is the form shown in Equation 115 and expected by pdepe
.
See Initial Value Problems for more information. For this
example, the resulting equation is
$${\pi}^{2}\frac{\partial u}{\partial t}={x}^{0}\frac{\partial}{\partial x}\left({x}^{0}\frac{\partial u}{\partial x}\right)+0$$
with parameter m = 0 and the terms
$$\begin{array}{l}c\left(x,t,u,\frac{\partial u}{\partial x}\right)={\pi}^{2}\\ f\left(x,t,u,\frac{\partial u}{\partial x}\right)=\frac{\partial u}{\partial x}\\ s\left(x,t,u,\frac{\partial u}{\partial x}\right)=0.\end{array}$$
Code the PDE. Once
you rewrite the PDE in the form shown above (Equation 115) and identify the terms, you
can code the PDE in a function that pdepe
can use.
The function must be of the form
[c,f,s] = pdefun(x,t,u,dudx)
where c
, f
, and s
correspond
to the c, f, and s terms.
The code below computes c
, f
,
and s
for the example problem.
function [c,f,s] = pdex1pde(x,t,u,DuDx) c = pi^2; f = DuDx; s = 0;
Code the initial conditions function. You must code the initial conditions in a function of the form
u = icfun(x)
The code below represents the initial conditions in the function pdex1ic
.
function u0 = pdex1ic(x) u0 = sin(pi*x);
Code the boundary conditions function. You must also code the boundary conditions in a function of the form
[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
The boundary conditions, written in the same form as Equation 117, are
$$u(0,t)+0\cdot \frac{\partial u}{\partial x}(0,t)=0\text{at}x=0$$
and
$$\pi {e}^{t}+1\cdot \frac{\partial u}{\partial x}(1,t)=0\text{at}x=1.$$
The code below evaluates the components p(x,t,u)
and q(x,t)
of the boundary conditions in the function pdex1bc
.
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = pi * exp(t); qr = 1;
In the function pdex1bc
, pl
and ql
correspond
to the left boundary conditions (x = 0), and pr
and qr
correspond
to the right boundary condition (x = 1).
Select mesh points
for the solution. Before you use the MATLAB PDE solver,
you need to specify the mesh points (t,x)
at which you want pdepe
to evaluate the solution.
Specify the points as vectors t
and x
.
The vectors t
and x
play
different roles in the solver (see PDE Solver).
In particular, the cost and the accuracy of the solution depend strongly
on the length of the vector x
. However, the computation
is much less sensitive to the values in the vector t
.
This example requests the solution on the mesh produced by 20
equally spaced points from the spatial interval [0,1] and five values
of t
from the time interval [0,2].
x = linspace(0,1,20); t = linspace(0,2,5);
Apply the PDE solver. The
example calls pdepe
with m = 0
,
the functions pdex1pde
, pdex1ic
,
and pdex1bc
, and the mesh defined by x
and t
at
which pdepe
is to evaluate the solution. The pdepe
function
returns the numerical solution in a threedimensional array sol
,
where sol(i,j,k)
approximates the k
th
component of the solution, u_{k},
evaluated at t(i)
and x(j)
.
m = 0; sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
This example uses @
to pass pdex1pde
, pdex1ic
,
and pdex1bc
as function handles to pdepe
.
Note: See Function Handles for more information about function handles. 
View the results. Complete the example by displaying the results:
Extract and display the first solution component. In this example, the solution u has only one component, but for illustrative purposes, the example "extracts" it from the threedimensional array. The surface plot shows the behavior of the solution.
u = sol(:,:,1); surf(x,t,u) title('Numerical solution computed with 20 mesh points') xlabel('Distance x') ylabel('Time t')
Display a solution profile at t_{f} , the final value of t. In this example, t_{f} = t = 2.
figure plot(x,u(end,:)) title('Solution at t = 2') xlabel('Distance x') ylabel('u(x,2)')
Evaluating the Solution. After obtaining and plotting the solution above, you might be
interested in a solution profile for a particular value of t
,
or the time changes of the solution at a particular point x
.
The k
th column u(:,k)
(of the
solution extracted in step 7) contains the time history
of the solution at x(k)
. The j
th
row u(j,:)
contains the solution profile at t(j)
.
Using the vectors x
and u(j,:)
,
and the helper function pdeval
,
you can evaluate the solution u
and its derivative
∂u/∂x at any set
of points xout
[uout,DuoutDx] = pdeval(m,x,u(j,:),xout)
The example pdex3
uses pdeval
to
evaluate the derivative of the solution at xout = 0
.
See pdeval
for details.
This example illustrates the solution of a system of partial differential equations. The problem is taken from electrodynamics. It has boundary layers at both ends of the interval, and the solution changes rapidly for small t.
The PDEs are
$$\begin{array}{l}\frac{\partial {u}_{1}}{\partial t}=0.024\frac{{\partial}^{2}{u}_{1}}{\partial {x}^{2}}F\left({u}_{1}{u}_{2}\right)\\ \frac{\partial {u}_{2}}{\partial t}=0.170\frac{{\partial}^{2}{u}_{2}}{\partial {x}^{2}}+F\left({u}_{1}{u}_{2}\right)\end{array}$$
where F(y) = exp(5.73y) – exp(–11.46y). The equations hold on an interval 0 ≤ x ≤ 1 for times t ≥ 0.
The solution u satisfies the initial conditions
$$\begin{array}{l}{u}_{1}(x,0)\equiv 1\\ {u}_{2}(x,0)\equiv 0\end{array}$$
and boundary conditions
$$\begin{array}{l}\frac{\partial {u}_{1}}{\partial x}(0,t)\equiv 0\\ {u}_{2}(0,t)\equiv 0\\ {u}_{1}(1,t)\equiv 1\\ \frac{\partial {u}_{2}}{\partial x}(1,t)\equiv 0.\end{array}$$
Note:
The file, 
Rewrite the PDE. In
the form expected by pdepe
, the equations are
$$\left[\begin{array}{c}1\\ 1\end{array}\right]\cdot \ast \frac{\partial}{\partial t}\left[\begin{array}{c}{u}_{1}\\ {u}_{2}\end{array}\right]=\frac{\partial}{\partial x}\left[\begin{array}{c}0.024\left(\partial {u}_{1}/\partial x\right)\\ 0.170\left(\partial {u}_{2}/\partial x\right)\end{array}\right]+\left[\begin{array}{c}F\left({u}_{1}{u}_{2}\right)\\ \text{}F\left({u}_{1}{u}_{2}\right)\end{array}\right].$$
The boundary conditions on the partial derivatives of u have
to be written in terms of the flux. In the form expected by pdepe
,
the left boundary condition is
$$\left[\begin{array}{c}0\\ {u}_{2}\end{array}\right]+\left[\begin{array}{c}1\\ 0\end{array}\right]\cdot \ast \left[\begin{array}{c}0.024\left(\partial {u}_{1}/\partial x\right)\\ 0.170\left(\partial {u}_{2}/\partial x\right)\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right]$$
and the right boundary condition is
$$\left[\begin{array}{c}{u}_{1}1\\ 0\end{array}\right]+\left[\begin{array}{c}0\\ 1\end{array}\right]\cdot \ast \left[\begin{array}{c}0.024\left(\partial {u}_{1}/\partial x\right)\\ 0.170\left(\partial {u}_{2}/\partial x\right)\end{array}\right]=\left[\begin{array}{c}0\\ 0\end{array}\right].$$
Code the PDE. After
you rewrite the PDE in the form shown above, you can code it as a
function that pdepe
can use. The function must
be of the form
[c,f,s] = pdefun(x,t,u,dudx)
where c
, f
, and s
correspond
to the c, f ,
and s terms in Equation 115.
function [c,f,s] = pdex4pde(x,t,u,DuDx) c = [1; 1]; f = [0.024; 0.17] .* DuDx; y = u(1)  u(2); F = exp(5.73*y)exp(11.47*y); s = [F; F];
Code the initial conditions function. The initial conditions function must be of the form
u = icfun(x)
The code below represents the initial conditions in the function pdex4ic
.
function u0 = pdex4ic(x); u0 = [1; 0];
Code the boundary conditions function. The boundary conditions functions must be of the form
[pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)
The code below evaluates the components p(x,t,u) and q(x,t) (Equation 117) of the boundary conditions in the
function pdex4bc
.
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t) pl = [0; ul(2)]; ql = [1; 0]; pr = [ur(1)1; 0]; qr = [0; 1];
Select mesh points for the solution. The solution changes rapidly for small t. The program selects the step size in time to resolve this sharp change, but to see this behavior in the plots, output times must be selected accordingly. There are boundary layers in the solution at both ends of [0,1], so mesh points must be placed there to resolve these sharp changes. Often some experimentation is needed to select the mesh that reveals the behavior of the solution.
x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1]; t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2];
Apply the PDE solver. The
example calls pdepe
with m = 0
,
the functions pdex4pde
, pdex4ic
,
and pdex4bc
, and the mesh defined by x and t
at which pdepe
is to evaluate the solution. The pdepe
function
returns the numerical solution in a threedimensional array sol
,
where sol(i,j,k)
approximates the k
th
component of the solution, μ_{k},
evaluated at t(i)
and x(j)
.
m = 0; sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
View the results. The surface plots show the behavior of the solution components.
u1 = sol(:,:,1); u2 = sol(:,:,2); figure surf(x,t,u1) title('u1(x,t)') xlabel('Distance x') ylabel('Time t')
figure surf(x,t,u2) title('u2(x,t)') xlabel('Distance x') ylabel('Time t')
The following additional examples are available. Type
edit examplename
to view the code and
examplename
to run the example.
Example Name  Description 

 Simple PDE that illustrates the straightforward formulation, computation, and plotting of the solution 
 Problem that involves discontinuities 
 Problem that requires computing values of the partial derivative 
 System of two PDEs whose solution has boundary layers at both ends of the interval and changes rapidly for small t. 
 System of PDEs with step functions as initial conditions 