# pdepe

Solve 1-D parabolic and elliptic PDEs

## Syntax

## Description

solves a system of parabolic and elliptic PDEs with one spatial variable
`sol`

= pdepe(`m`

,`pdefun`

,`icfun`

,`bcfun`

,`xmesh`

,`tspan`

)*x* and time *t*. At least one equation must be
parabolic. The scalar `m`

represents the symmetry of the problem (slab,
cylindrical, or spherical). The equations being solved are coded in
`pdefun`

, the initial value is coded in `icfun`

, and the
boundary conditions are coded in `bcfun`

. The ordinary differential
equations (ODEs) resulting from discretization in space are integrated to obtain approximate
solutions at the times specified in `tspan`

. The `pdepe`

function returns values of the solution on a mesh provided in
`xmesh`

.

`[`

also finds where functions of (`sol`

,`tsol`

,`sole`

,`te`

,`ie`

] = pdepe(`m`

,`pdefun`

,`icfun`

,`bcfun`

,`xmesh`

,`tspan`

,`options`

)*t*,*u*(*x*,*t*)), called event functions, are zero. In the output, `te`

is
the time of the event, `sole`

is the solution at the time of the event, and
`ie`

is the index of the triggered event. `tsol`

is a
column vector of times specified in `tspan`

, prior to the first terminal
event.

For each event function, specify whether the integration is to terminate at a zero and
whether the direction of the zero-crossing matters. Do this by setting the
`'Events'`

option of `odeset`

to a function, such as
`@myEventFcn`

, and creating a corresponding function:
[`value`

,`isterminal`

,`direction`

] =
`myEventFcn`

(`m`

,`t`

,`xmesh`

,`umesh`

).
The `xmesh`

input contains the spatial mesh and `umesh`

is
the solution at the mesh points.

## Examples

### Solve Heat Equation in Cylindrical Coordinates

Solve the heat equation in cylindrical coordinates using `pdepe`

, and plot the solution.

In cylindrical coordinates with angular symmetry the heat equation is

$$\frac{\partial \mathit{u}}{\partial \mathit{t}}=\frac{1}{\mathit{x}}\frac{\partial}{\partial \mathit{x}}\left(\mathit{x}\frac{\partial \mathit{u}}{\partial \mathit{x}}\right).$$

The equation is defined for $0\le \mathit{x}\le 1$ at times $\mathit{t}\ge 0$. The initial condition is defined in terms of the bessel function ${\mathit{J}}_{0}\left(\mathit{x}\right)$ and its first zero $\mathit{n}=2.404825557695773$ as

$$\mathit{u}\left(\mathit{x},0\right)={\mathit{J}}_{0}\left(\mathrm{nx}\right).$$

Since this problem is in cylindrical coordinates (`m = 1`

), `pdepe`

automatically enforces the symmetry condition at $\mathit{x}=0$. The right boundary condition is

$$\mathit{u}\left(1,\mathit{t}\right)={\mathit{J}}_{0}\left(\mathit{n}\right){\mathit{e}}^{-{\mathit{n}}^{2}\mathit{t}}.$$

The initial and boundary conditions are selected to be consistent with the analytic solution to the problem,

$$\mathit{u}\left(\mathit{x},\mathit{t}\right)={\mathit{J}}_{0}\left(\mathrm{nx}\right){\mathit{e}}^{-{\mathit{n}}^{2}\mathit{t}}.$$

To solve this equation in MATLAB®, you need to code the equation, initial conditions, and boundary conditions, then select a suitable solution mesh before calling the solver `pdepe`

. You either can include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.

**Code Equation**

Before you can code the equation, you need to rewrite it in a form that the `pdepe`

solver expects. The standard form that `pdepe`

expects is

$$\mathit{c}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)\frac{\partial \mathit{u}}{\partial \mathit{t}}={\mathit{x}}^{-\mathit{m}}\frac{\partial}{\partial \mathit{x}}\left({\mathit{x}}^{\mathit{m}}\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)\right)+\mathit{s}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right).$$

Written in this form, the PDE becomes

$$\frac{\partial \mathit{u}}{\partial \mathit{t}}={\mathit{x}}^{-1}\frac{\partial}{\partial \mathit{x}}\left({\mathit{x}}^{1}\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)$$

With the equation in the proper form you can read off the relevant terms:

$\mathit{m}=1$

$\mathit{c}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=1$

$\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=\frac{\partial \mathit{u}}{\partial \mathit{x}}$

$\mathit{s}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0$

Now you can create a function to code the equation. The function should have the signature `[c,f,s] = heatcyl(x,t,u,dudx)`

:

`x`

is the independent spatial variable.`t`

is the independent time variable.`u`

is the dependent variable being differentiated with respect to`x`

and`t`

.`dudx`

is the partial spatial derivative $\partial \mathit{u}/\partial \mathit{x}$.The outputs

`c`

,`f`

, and`s`

correspond to coefficients in the standard PDE equation form expected by`pdepe`

. These coefficients are coded in terms of the input variables`x`

,`t`

,`u`

, and`dudx`

.

As a result, the equation in this example can be represented by the function:

function [c,f,s] = heatcyl(x,t,u,dudx) c = 1; f = dudx; s = 0; end

(Note: All functions are included as local functions at the end of the example.)

**Code Initial Conditions**

Next, write a function that returns the initial condition. The initial condition is applied at the first time value `tspan(1)`

. The function should have the signature `u0 = heatic(x)`

.

The corresponding function for $\mathit{u}\left(\mathit{x},0\right)={\mathit{J}}_{0}\left(\mathrm{nx}\right)$ is

function u0 = heatic(x) n = 2.404825557695773; u0 = besselj(0,n*x); end

**Code Boundary Conditions**

Now, write a function that evaluates the boundary condition

$\mathit{u}\left(1,\mathit{t}\right)={\mathit{J}}_{0}\left(\mathit{n}\right){\mathit{e}}^{-{\mathit{n}}^{2}\mathit{t}}.$

Since this problem is in cylindrical coordinates (`m = 1`

), `pdepe`

automatically enforces the symmetry condition at $\mathit{x}=0$, so you do not need to specify a left boundary condition.

For problems posed on the interval $\mathit{a}\le \mathit{x}\le \mathit{b}$, the boundary conditions apply for all $\mathit{t}$ and either $\mathit{x}=\mathit{a}$ or $\mathit{x}=\mathit{b}$. The standard form for the boundary conditions expected by the solver is

$$\mathit{p}\left(\mathit{x},\mathit{t},\mathit{u}\right)+\mathit{q}\left(\mathit{x},\mathit{t}\right)\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0.$$

Written in this form, the boundary conditions for the partial derivatives of $\mathit{u}$ need to be expressed in terms of the flux $\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)$. So the right boundary condition for this problem is

$$\mathit{u}\left(1,\mathit{t}\right)-{\mathit{J}}_{0}\left(\mathit{n}\right){\mathit{e}}^{-{\mathit{n}}^{2}\mathit{t}}+0\cdot \mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0.$$

The boundary function should have the function signature `[pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)`

:

The inputs

`xl`

and`ul`

correspond to $\mathit{u}$ and $\mathit{x}$ for the left boundary.The inputs

`xr`

and`ur`

correspond to $\mathit{u}$ and $\mathit{x}$ for the right boundary.`t`

is the independent time variable.The outputs

`pl`

and`ql`

correspond to ${\mathit{p}}_{\mathit{L}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$ and ${\mathit{q}}_{\mathit{L}}\left(\mathit{x},\mathit{t}\right)$ for the left boundary ($\mathit{x}=0$ for this problem).The outputs

`pr`

and`qr`

correspond to ${\mathit{p}}_{\mathit{R}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$ and ${\mathit{q}}_{\mathit{R}}\left(\mathit{x},\mathit{t}\right)$ for the right boundary ($\mathit{x}=1$ for this problem).

The boundary conditions in this example are represented by the function:

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t) n = 2.404825557695773; pl = 0; %ignored by solver since m=1 ql = 0; %ignored by solver since m=1 pr = ur-besselj(0,n)*exp(-n^2*t); qr = 0; end

**Select Solution Mesh**

Before solving the equation you need to specify the mesh points $\left(\mathit{t},\mathit{x}\right)$ 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. In particular, the cost and 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`

.

For this problem, use a mesh with 25 equally spaced points in the spatial interval [0,1] for both `x`

and `t`

.

x = linspace(0,1,25); t = linspace(0,1,25);

**Solve Equation**

Finally, solve the equation using the symmetry `m`

, the PDE equation, the initial conditions, the boundary conditions, and the meshes for `x`

and `t`

.

m = 1; sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,t);

`pdepe`

returns the solution in a 3-D array `sol`

, where `sol(i,j,k)`

approximates the `k`

th component of the solution ${\mathit{u}}_{\mathit{k}}$ evaluated at `t(i)`

and `x(j)`

. The size of `sol`

is `length(t)`

-by-`length(x)`

-by-`length(u0)`

, since `u0`

specifies an initial condition for each solution component. For this problem, `u`

has only one component, so `sol`

is a 25-by-25 matrix, but in general you can extract the `k`

th solution component with the command `u = sol(:,:,k)`

.

Extract the first solution component from `sol`

.

u = sol(:,:,1);

**Plot Solution**

Create a surface plot of the solution. Since the problem is posed with cylindrical coordinates on a disc, the `x`

values show the temperature on the disc some distance from the center, and the `t`

values show how the temperature at a particular location changes over time.

surf(x,t,u) xlabel('x') ylabel('t') zlabel('u(x,t)') view([150 25])

Plot the temperature change at the center ($\mathit{x}=0$) of the disc.

plot(t,sol(:,1)) xlabel('Time') ylabel('Temperature u(0,t)') title('Temperature change at center of disc')

**Local Functions**

Listed here are the local helper functions that the PDE solver `pdepe`

calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

function [c,f,s] = heatcyl(x,t,u,dudx) c = 1; f = dudx; s = 0; end %---------------------------------------------- function u0 = heatic(x) n = 2.404825557695773; u0 = besselj(0,n*x); end %---------------------------------------------- function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t) n = 2.404825557695773; pl = 0; %ignored by solver since m=1 ql = 0; %ignored by solver since m=1 pr = ur-besselj(0,n)*exp(-n^2*t); qr = 0; end %----------------------------------------------

### Solve Oscillatory PDE with Event Logging

Solve a partial differential equation and use an event function to log zero-crossings in the oscillatory solution.

Consider the equation

$$\frac{1}{\mathit{x}}\frac{\partial \mathit{u}}{\partial \mathit{t}}=\frac{\partial}{\partial \mathit{x}}\left(\frac{1}{\mathit{t}}\mathit{u}\right).$$

The equation is defined for $0\le \mathit{x}\le 1$ and $\mathit{t}\ge \mathrm{0}.1$. The initial condition is

$$\mathit{u}\left(\mathit{x},0.1\right)=1.$$

The boundary conditions are

$$\mathit{u}\left(0,\mathit{t}\right)=1,$$

$$\mathit{u}\left(1,\mathit{t}\right)=\mathrm{cos}\left(\pi \text{\hspace{0.17em}}\mathit{t}\right).$$

Additionally, the zero-crossings of the solution are of interest.

To solve this equation in MATLAB, you need to code the equation, initial conditions, boundary conditions, and event function, then select a suitable solution mesh before calling the solver `pdepe`

. You either can include the required functions as local functions at the end of a file (as done here), or save them as separate, named files in a directory on the MATLAB path.

**Code Equation**

Before you can code the equation, you need to rewrite it in a form that the `pdepe`

solver expects. The standard form that `pdepe`

expects is

$$\mathit{c}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)\frac{\partial \mathit{u}}{\partial \mathit{t}}={\mathit{x}}^{-\mathit{m}}\frac{\partial}{\partial \mathit{x}}\left({\mathit{x}}^{\mathit{m}}\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)\right)+\mathit{s}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right).$$

The PDE equation is already in this form:

$\frac{1}{\mathit{x}}\frac{\partial \mathit{u}}{\partial \mathit{t}}=\frac{\partial}{\partial \mathit{x}}\left(\frac{1}{\mathit{t}}\mathit{u}\right).$

So you can read off the relevant terms:

$$\mathit{m}=0$$

$$\mathit{c}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=\frac{1}{\mathit{x}}$$

$$\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=\frac{1}{\mathit{t}}\mathit{u}$$

$$\mathit{s}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0$$

Now you can create a function to code the equation. The function should have the signature `[c,f,s] = oscpde(x,t,u,dudx)`

:

`x`

is the independent spatial variable.`t`

is the independent time variable.`u`

is the dependent variable being differentiated with respect to`x`

and`t`

.`dudx`

is the partial spatial derivative $\partial \mathit{u}/\partial \mathit{x}$.The outputs

`c`

,`f`

, and`s`

correspond to coefficients in the standard PDE equation form expected by`pdepe`

. These coefficients are coded in terms of the input variables`x`

,`t`

,`u`

, and`dudx`

.

As a result, the equation in this example can be represented by the function:

function [c,f,s] = oscpde(x,t,u,dudx) c = 1/x; f = u/t; s = 0; end

(Note: All functions are included as local functions at the end of the example.)

**Code Initial Conditions**

Next, write a function that returns the initial condition. The initial condition is applied at the first time value `tspan(1)`

. The function should have the signature `u0 = oscic(x)`

.

The corresponding function for $\mathit{u}\left(\mathit{x},0.1\right)=1$ is

function u0 = oscic(x) u0 = 1; end

**Code Boundary Conditions**

Now, write a function that evaluates the boundary conditions

$$\mathit{u}\left(0,\mathit{t}\right)=1,$$

$$\mathit{u}\left(1,\mathit{t}\right)=\mathrm{cos}\left(\pi \text{\hspace{0.17em}}\mathit{t}\right).$$

For problems posed on the interval $\mathit{a}\le \mathit{x}\le \mathit{b}$, the boundary conditions apply for all $\mathit{t}$ and either $\mathit{x}=\mathit{a}$ or $\mathit{x}=\mathit{b}$. The standard form for the boundary conditions expected by the solver is

$$\mathit{p}\left(\mathit{x},\mathit{t},\mathit{u}\right)+\mathit{q}\left(\mathit{x},\mathit{t}\right)\mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0.$$

Written in this form, the boundary conditions for this problem become

$$\mathit{u}\left(0,\mathit{t}\right)-1+0\cdot \mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0,$$

$$\mathit{u}\left(1,\mathit{t}\right)-\mathrm{cos}\left(\pi \text{\hspace{0.17em}}\mathit{t}\right)+0\cdot \mathit{f}\left(\mathit{x},\mathit{t},\mathit{u},\frac{\partial \mathit{u}}{\partial \mathit{x}}\right)=0.$$

The boundary function should have the function signature `[pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t)`

:

The inputs

`xl`

and`ul`

correspond to $\mathit{x}$ and $\mathit{u}$ for the left boundary.The inputs

`xr`

and`ur`

correspond to $\mathit{x}$ and $\mathit{u}$ for the right boundary.`t`

is the independent time variable.The outputs

`pl`

and`ql`

correspond to ${\mathit{p}}_{\mathit{L}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$ and ${\mathit{q}}_{\mathit{L}}\left(\mathit{x},\mathit{t}\right)$ for the left boundary ($\mathit{x}=0$ for this problem).The outputs

`pr`

and`qr`

correspond to ${\mathit{p}}_{\mathit{R}}\left(\mathit{x},\mathit{t},\mathit{u}\right)$ and ${\mathit{q}}_{\mathit{R}}\left(\mathit{x},\mathit{t}\right)$ for the right boundary ($\mathit{x}=1$ for this problem).

The boundary conditions in this example are represented by the function:

function [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t) pl = ul - 1; ql = 0; pr = ur - cos(pi*t); qr = 0; end

**Code Event Function**

Use an event function to log the zero-crossings of the solution in the integration. The event function has a function signature of `[value,isterminal,direction] = pdevents(m,t,xmesh,umesh)`

:

`m`

is the coordinate symmetry specified as the first input to`pdepe`

.`t`

is the current time (a scalar).`xmesh`

is the spatial mesh.`umesh`

contains the solution at the mesh points.`value`

is an equation of interest, usually expressed in terms of the solution`umesh`

. When`value`

equals 0, an event occurs.`isterminal`

specifies whether an event stops the integration. If`isterminal`

is 0, events are logged but the integration does not stop. If`isterminal`

is 1, then the integration stops when an event occurs.`direction`

specifies the direction of zero crossing. If 1, only zero-crossings with a positive slope trigger an event. If -1, the zero-crossings must have a negative slope. If 0, then any zero-crossing triggers an event.

At each time in the integration, the solver calls the event function to check for zero crossings. To log all zero crossings, `value`

should look for changes in sign in the solution vector `umesh`

. Specify `isterminal`

and `direction`

as vectors of zeros with the same size, since the events in this example are not terminal and the zero-crossings can occur with any slope.

The event function for this problem is

function [value,isterminal,direction] = pdevents(m,t,xmesh,umesh) value = umesh; isterminal = zeros(size(umesh)); direction = zeros(size(umesh)); end

**Select Solution Mesh**

Before solving the equation you need to specify the mesh points $\left(\mathit{t},\mathit{x}\right)$ at which you want `pdepe`

to evaluate the solution. For this problem, use a fine mesh of 50 points in the intervals $0\le \mathit{x}\le 1$ and $0.1\le \mathit{t}\le 1$. A fine mesh gives good resolution of the oscillatory solution.

x = linspace(0,1,50); t = linspace(0.1,pi,50);

**Solve Equation**

Finally, solve the equation using the symmetry `m`

, the PDE equation, the initial conditions, the boundary conditions, the event function, and the meshes for `x`

and `t`

. Use `odeset`

to create an options structure that references the events function, and pass in the structure as the last input argument to `pdepe`

. Specify five output arguments to return information from both the event function and the solver:

`sol`

is the solution computed by`pdepe`

.`tsol`

is a vector of times before a terminal event.`tsol`

is equal to`t`

when no events are terminal.`sole`

is the solution at the time of each event.`te`

is the time of each event.`ie`

is the index of each event. Since`values = umesh`

in the event function,`ie`

gives the indices of`umesh`

that triggered an event at each time step.

```
m = 0;
options = odeset('Events',@pdevents);
[sol,tsol,sole,te,ie] = pdepe(m,@oscpde,@oscic,@oscbc,x,t,options);
```

Extract the solution as a matrix `u`

.

u = sol(:,:,1);

**Plot Solution**

Create a surface plot of the solution and view the plot from above.

surf(x,t,u) view(2)

Plot the points where events occurred, with a surface $\mathit{f}\left(\mathit{x},\mathit{t}\right)=0$ for reference. The output index vector `ie`

is useful to pick out the event locations. The expression `x(ie)'`

gives the x-values where events occurred, and the expression `sole(x==x(ie)')`

gives the corresponding solution values.

view([39 30]) xlabel('x') ylabel('t') zlabel('u(x,t)') hold on plot3(x(ie)',te,sole(x==x(ie)'),'r*') surf(x,t,zeros(size(u)),'EdgeColor','flat') hold off

**Local Functions**

Listed here are the local helper functions that the PDE solver pdepe calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

function [c,f,s] = oscpde(x,t,u,dudx) c = 1/x; f = u/t; s = 0; end %---------------------------------------------- function u0 = oscic(x) u0 = 1; end %---------------------------------------------- function [pl,ql,pr,qr] = oscbc(xl,ul,xr,ur,t) pl = ul - 1; ql = 0; pr = ur - cos(pi*t); qr = 0; end %---------------------------------------------- function [value, isterminal, direction] = pdevents(m,t,xmesh,umesh) value = umesh; isterminal = zeros(size(umesh)); direction = zeros(size(umesh)); end %----------------------------------------------

## Input Arguments

`m`

— Symmetry constant

`0`

| `1`

| `2`

Symmetry constant, specified as one of the values in this table.
`m`

specifies the type of problem in the equations specified by
`pdefun`

. Once you rewrite the equations in the form expected by
`pdepe`

, you can read the value of `m`

from the
equation.

Value | Description | Laplacian of flux |
---|---|---|

| 1-D Cartesian coordinates with no symmetry | $$\Delta f=\frac{{\partial}^{2}f}{\partial {x}^{2}}$$ |

| 2-D Cylindrical coordinates with azimuthal angular symmetry | $$\Delta f=\frac{1}{\rho}\frac{\partial}{\partial \rho}\left(\rho \frac{\partial f}{\partial \rho}\right)+\frac{1}{{\rho}^{2}}\frac{{\partial}^{2}f}{\partial {\phi}^{2}}$$, where $$\frac{\partial f}{\partial \phi}=0$$ (azimuthal symmetry). |

| 3-D Spherical coordinates with azimuthal and zenith angular symmetries | $$\Delta f=\frac{1}{{r}^{2}}\frac{\partial}{\partial r}\left({r}^{2}\frac{\partial f}{\partial r}\right)+\frac{1}{{r}^{2}\mathrm{sin}\theta}\frac{\partial}{\partial \theta}\left(\mathrm{sin}\theta \frac{\partial f}{\partial \theta}\right)+\frac{1}{{r}^{2}{\mathrm{sin}}^{2}\theta}\frac{{\partial}^{2}f}{\partial {\phi}^{2}}$$, where $$\frac{\partial f}{\partial \theta}=\frac{\partial f}{\partial \phi}=0$$ (azimuthal and zenith angular symmetries). |

When `m > 0`

, `pdepe`

requires the
left spatial boundary *a* to be ≥ 0, and it imposes the left boundary
condition automatically to account for the coordinate discontinuity at the origin. In
that case, `pdepe`

ignores any specified left boundary conditions in
`bcfun`

.

`pdefun`

— PDE function for equations

function handle

PDE function for equations, specified as a function handle that defines the coefficients of the PDE.

`pdefun`

encodes the coefficients for PDEs 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).$$

The terms in the equation are:

$$f\left(x,t,u,\frac{\partial u}{\partial x}\right)$$ is a flux term.

$$s\left(x,t,u,\frac{\partial u}{\partial x}\right)$$ is a source term.

The coupling of the partial derivatives with respect to time is restricted to multiplication by a diagonal matrix $$c\left(x,t,u,\frac{\partial u}{\partial x}\right)$$. The diagonal elements of this matrix are either zero or positive. An element that is zero corresponds to an elliptic equation, and all other elements correspond to parabolic equations. 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 those values of*x*are mesh points.Discontinuities in

*c*and*s*due to material interfaces are permitted provided that a mesh point is placed at each interface.

The PDEs hold for *t*_{0} ≤ *t* ≤
*t _{f}* and

*a*≤

*x*≤

*b*. The values

`tspan(1)`

and
`tspan(end)`

correspond to
*t*

_{0}and

*t*, while

_{f}`xmesh(1)`

and
`xmesh(end)`

correspond to *a*and

*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*must be ≥ 0.

#### Coding the Equation

Write a function `pdefun`

that calculates the values of the
coefficients *c*, *f*, and *s*. Use
the function
signature

[c,f,s] = pdefun(x,t,u,dudx)

The input arguments to `pdefun`

are scalars `x`

and `t`

and vectors `u`

and `dudx`

.
The vector `u`

approximates the solution *u*, and
`dudx`

approximates its partial derivative with respect to
*x*. The output arguments `c`

,
`f`

, and `s`

are column vectors with a number of
elements equal to the number of equations. `c`

stores the diagonal
elements of the matrix *c*.

#### Example

The heat equation $$\frac{\partial u}{\partial t}=\frac{{\partial}^{2}u}{\partial {x}^{2}}$$ can be rewritten in the form expected by the solver as

$$1\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\frac{\partial u}{\partial t}={x}^{0}\frac{\partial}{\partial x}\left({x}^{0}\frac{\partial u}{\partial x}\right).$$

In this form you can read off the value of the coefficients as *m* = 0, *c* = 1, *f* =
∂*u*/∂*x*, and *s* = 0. The function for this equation is:

function [c,f,s] = heatpde(x,t,u,dudx) c = 1; f = dudx; s = 0; end

**Data Types: **`function_handle`

`icfun`

— Initial condition function

function handle

Initial condition function, specified as a function handle that defines the initial condition.

For *t* = *t*_{0} and all
*x*, the solution components satisfy initial conditions of the form

$$u\left(x,{t}_{0}\right)={u}_{0}\left(x\right).$$

#### Coding the Initial Condition

Write a function `icfun`

that defines the initial conditions. Use
the function signature:

function u0 = icfun(x)

When called with an argument `x`

, the function
`icfun`

evaluates and returns the initial values of the solution
components at `x`

in the column vector `u0`

. The
number of elements in `u0`

is equal to the number of
equations.

#### Example

The constant initial condition $$u\left(x,0\right)=1/2$$ corresponds to this function:

function u0 = heatic(x) u0 = 0.5; end

**Data Types: **`function_handle`

`bcfun`

— Boundary condition function

function handle

Boundary condition function, specified as a function handle that defines the boundary conditions.

For all *t* and either *x* = *a*
or *x* = *b*, the solution components satisfy a
boundary condition of the form

$$p\left(x,t,u\right)+q\left(x,t\right)f\left(x,t,u,\frac{\partial u}{\partial x}\right)=0.$$

Elements of *q* are either zero or never zero. Note that the
boundary conditions are expressed in terms of the flux *f* rather than
∂*u*/∂*x*. Also, of the two coefficients, only
*p* can depend on *u*.

#### Coding the Boundary Conditions

Write a function `bcfun`

that defines the terms
*p* and *q* of the boundary conditions. Use the
function signature

function [pl,ql,pr,qr] = bcfun(xl,ul,xr,ur,t)

`ul`

is the approximate solution at the left boundary`xl = a`

, and`ur`

is the approximate solution at the right boundary`xr = b`

.`pl`

and`ql`

are column vectors corresponding to*p*and*q*evaluated at`xl`

. Similarly,`pr`

and`qr`

correspond to`xr`

.The number of elements in the outputs

`pl`

,`ql`

,`pr`

, and`qr`

is equal to the number of equations.

When *m* > 0 and *a* = 0, boundedness of the solution near *x* = 0 requires that the flux *f* vanish at *a* = 0. `pdepe`

imposes this boundary condition
automatically, and it ignores values returned in `pl`

and
`ql`

.

#### Example

Consider the boundary conditions $$u\left(0,t\right)=0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}u\left(1,t\right)=1$$ on the interval $$0\le x\le 1$$. In the form expected by the solver, the boundary conditions are

$$\begin{array}{l}u\left(0,t\right)+0\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}f\left(0,t,u,\frac{\partial u}{\partial x}\right)=0,\\ u\left(1,t\right)-1+0\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}f\left(0,t,u,\frac{\partial u}{\partial x}\right)=0.\end{array}$$

In this form it is clear that both *q* terms are zero. The
*p* terms are nonzero to reflect the values of
*u*. The corresponding function is

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur - 1; qr = 0; end

**Data Types: **`function_handle`

`xmesh`

— Spatial mesh

vector

Spatial mesh, specified as a vector `[x0 x1 ... xn]`

specifying the
points at which a numerical solution is requested for every value in
`tspan`

. Second-order approximations to the solution are made on the
mesh specified in `xmesh`

. `pdepe`

does
*not* select the mesh in *x* automatically. You
must provide an appropriate fixed mesh in `xmesh`

.

The elements of `xmesh`

must satisfy ```
x0 < x1 < ...
< xn
```

. Generally, it is best to use closely spaced mesh points where the
solution changes rapidly. The length of `xmesh`

must be
`>=3`

, and the computational cost of the solution depends strongly on
the length of `xmesh`

.

For problems with discontinuities, you should place mesh points at the
discontinuities so that the problem is smooth within each subinterval. However, when *m* > 0 , it is not necessary to use a fine mesh near *x* = 0 to account for the coordinate singularity.

**Example: **`xmesh = linspace(0,5,25)`

uses a spatial mesh of 25 points
between 0 and 5.

**Data Types: **`single`

| `double`

`tspan`

— Time span of integration

vector

Time span of integration, specified as a vector `[t0 t1 ... tf]`

specifying the points at which a solution is requested for every value in
`xmesh`

. The `pdepe`

function performs the time
integration with an ODE solver that selects both the time step and the formula
dynamically. The elements of `tspan`

merely specify where you want
answers, and the computational cost of the solution therefore depends only weakly on the
length of `tspan`

.

The elements of `tspan`

must satisfy ```
t0 < t1 < ...
< tf
```

. The length of `tspan`

must be
`>=3`

.

**Example: **`tspan = linspace(0,5,5)`

uses five time points between 0
and 5.

**Data Types: **`single`

| `double`

`options`

— Option structure

structure array

Option structure, specified as a structure array. Use the `odeset`

function to create or modify the option structure.
`pdepe`

supports these options:

In most cases, default values for these options provide satisfactory solutions.

**Example: **`options = odeset('RelTol',1e-5)`

specifies a relative
error tolerance of `1e-5`

.

**Data Types: **`struct`

## Output Arguments

`sol`

— Solution array

3-D array

Solution array, returned as a 3-D array.

`pdepe`

returns the solution as a multidimensional array.
*u _{i}* =

`ui`

= `sol`

(`:`

,`:`

,`i`

)
is an approximation to the `i`

th component of the solution vector
*u*. The element

`ui`

(`j`

,`k`

) =
`sol`

(`j`

,`k`

,`i`

)
approximates *u*at (

_{i}*t*,

*x*) = (

`tspan`

(`j`

),`xmesh`

(`k`

)).`ui`

=
`sol`

(`j`

,`:`

,`i`

)
approximates component `i`

of the solution at time
`tspan`

(`j`

) and mesh points
`xmesh(:)`

. Use `pdeval`

to compute the
approximation and its partial derivative
∂*u _{i}*/∂

*x*at points not included in

`xmesh`

. See `pdeval`

for details.`tsol`

— Solution times

column vector

Solution times, returned as a column vector of times specified in
`tspan`

that are prior to the first terminal event.

`sole`

— Solution at time of events

array

Solution at time of events, returned as an array. The event times in
`te`

correspond to solutions returned in `sole`

, and
`ie`

specifies which event occurred.

`te`

— Time of events

column vector

Time of events, returned as a column vector. The event times in
`te`

correspond to solutions returned in `sole`

, and
`ie`

specifies which event occurred.

`ie`

— Index of triggered event function

column vector

Index of triggered event function, returned as a column vector. The event times in
`te`

correspond to solutions returned in `sole`

, and
`ie`

specifies which event occurred.

## Tips

If

`uji = sol(j,:,i)`

approximates component`i`

of the solution at time`tspan(j)`

and mesh points`xmesh`

, then`pdeval`

evaluates the approximation and its partial derivative ∂*u*/∂_{i}*x*at the array of points`xout`

and returns them in`uout`

and`duoutdx`

:`[uout,duoutdx] = pdeval(m,xmesh,uji,xout)`

. The`pdeval`

function evaluates the partial derivative ∂*u*/∂_{i}*x*rather than the flux. The flux is continuous, but at a material interface the partial derivative may have a jump.

## Algorithms

The time integration is done with the `ode15s`

solver. `pdepe`

exploits the capabilities of
`ode15s`

for solving the differential-algebraic equations that arise when
the PDE contains elliptic equations, and for handling Jacobians with a specified sparsity
pattern.

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.

## References

[1] Skeel, R. D. and M. Berzins, "A
Method for the Spatial Discretization of Parabolic Equations in One Space Variable,"
*SIAM Journal on Scientific and Statistical Computing*, Vol. 11, 1990,
pp.1–32.

## Extended Capabilities

### Thread-Based Environment

Run code in the background using MATLAB® `backgroundPool`

or accelerate code with Parallel Computing Toolbox™ `ThreadPool`

.

## Version History

**Introduced before R2006a**

## Open Example

You have a modified version of this example. Do you want to open this example with your edits?

## MATLAB Command

You clicked a link that corresponds to this MATLAB command:

Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)