Before you create boundary conditions, you need to create a `PDEModel`

container.
For details, see Solve Problems Using Legacy PDEModel Objects. Suppose
that you have a container named `model`

, and that
the geometry is stored in `model`

. Examine the geometry
to see the label of each edge or face.

pdegplot(model,'EdgeLabels','on') % for 2-D pdegplot(model,'FaceLabels','on') % for 3-D

Now you can specify the boundary conditions for each edge or
face. If you have a system of PDEs, you can set a different boundary
condition for each component on each boundary edge or face. If the
boundary condition is a function of position, time, or the solution *u*,
set boundary conditions by using the syntax in Nonconstant Boundary Conditions.

If you do not specify a boundary condition for an edge or face,
the default is the Neumann boundary condition with the zero values
for `'g'`

and `'q'`

.

The Dirichlet boundary condition implies that the solution *u* on
a particular edge or face satisfies the equation

*hu* = *r*,

where *h* and *r* are functions
defined on ∂Ω, and can be functions of space (*x*, *y*,
and, in 3-D, *z*), the solution *u*,
and, for time-dependent equations, time. Often, you take *h* = 1, and set *r* to
the appropriate value. You can specify Dirichlet boundary conditions
as the value of the solution `u`

on the boundary
or as a pair of the parameters `h`

and `r`

.

Suppose that you have a PDE model named `model`

,
and edges or faces `[e1,e2,e3]`

, where the solution *u* must
equal `2`

. Specify this boundary condition as follows.

% For 3-D geometry: applyBoundaryCondition(model,'dirichlet','Face',[e1,e2,e3],'u',2); % For 2-D geometry: applyBoundaryCondition(model,'dirichlet','Edge',[e1,e2,e3],'u',2);

If the solution on edges or faces `[e1,e2,e3]`

satisfies
the equation 2*u* = 3, specify the boundary condition
as follows.

% For 3-D geometry: applyBoundaryCondition(model,'dirichlet','Face',[e1,e2,e3],'r',3,'h',2); % For 2-D geometry: applyBoundaryCondition(model,'dirichlet','Edge',[e1,e2,e3],'r',3,'h',2);

If you do not specify

`'r'`

,`applyBoundaryCondition`

sets its value to`0`

.If you do not specify

`'h'`

,`applyBoundaryCondition`

sets its value to`1`

.

The Dirichlet boundary condition for a system of PDEs is **hu** = **r**, where **h** is
a matrix, **u** is the solution vector,
and **r** is a vector.

Suppose that you have a PDE model named `model`

,
and edge or face labels `[e1,e2,e3]`

where the first
component of the solution *u* must equal `1`

,
while the second and third components must equal `2`

.
Specify this boundary condition as follows.

% For 3-D geometry: applyBoundaryCondition(model,'dirichlet','Face',[e1,e2,e3],... 'u',[1,2,2],'EquationIndex',[1,2,3]); % For 2-D geometry: applyBoundaryCondition(model,'dirichlet','Edge',[e1,e2,e3],... 'u',[1,2,2],'EquationIndex',[1,2,3]);

The

`'u'`

and`'EquationIndex'`

arguments must have the same length.If you exclude the

`'EquationIndex'`

argument, the`'u'`

argument must have length*N*.If you exclude the

`'u'`

argument,`applyBoundaryCondition`

sets the components in`'EquationIndex'`

to`0`

.

Suppose that you have a PDE model named `model`

,
and edge or face labels `[e1,e2,e3]`

where the first,
second, and third components of the solution *u* must
satisfy the equations 2*u _{1}* = 3, 4

H0 = [2 0 0; 0 4 0; 0 0 6]; R0 = [3;5;7]; % For 3-D geometry: applyBoundaryCondition(model,'dirichlet', ... 'Face',[e1,e2,e3], ... 'h',H0,'r',R0); % For 2-D geometry: applyBoundaryCondition(model,'dirichlet', ... 'Edge',[e1,e2,e3], ... 'h',H0,'r',R0);

The

`'r'`

parameter must be a numeric vector of length*N*. If you do not specify`'r'`

,`applyBoundaryCondition`

sets the values to`0`

.The

`'h'`

parameter can be an*N*-by-*N*numeric matrix or a vector of length*N*^{2}corresponding to the linear indexing form of the*N*-by-*N*matrix. For details about the linear indexing form, see Array Indexing (MATLAB). If you do not specify`'h'`

,`applyBoundaryCondition`

sets the value to the identity matrix.

Generalized Neumann boundary conditions imply that the solution *u* on
the edge or face satisfies the equation

$$\overrightarrow{n}\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\nabla u\right)+qu=g$$

The coefficient *c* is the same as the coefficient
of the second-order differential operator in the PDE equation

$$-\nabla \cdot \left(c\nabla u\right)+au=f\text{ondomain}\Omega $$

$$\overrightarrow{n}$$ is the outward unit normal. *q* and *g* are
functions defined on ∂Ω, and can be functions of space
(*x*, *y*, and, in 3-D, *z*),
the solution *u*, and, for time-dependent equations,
time.

Suppose that you have a PDE model named `model`

,
and edges or faces `[e1,e2,e3]`

where the solution *u* must
satisfy the Neumann boundary condition with `q = 2`

and ```
g
= 3
```

. Specify this boundary condition as follows.

% For 3-D geometry: applyBoundaryCondition(model,'neumann','Face',[e1,e2,e3],'q',2,'g',3); % For 2-D geometry: applyBoundaryCondition(model,'neumann','Edge',[e1,e2,e3],'q',2,'g',3);

If you do not specify

`'g'`

,`applyBoundaryCondition`

sets its value to`0`

.If you do not specify

`'q'`

,`applyBoundaryCondition`

sets its value to`0`

.

Neumann boundary conditions for a system of PDEs is $$n\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)+qu=g$$. For 2-D systems, the notation $$n\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)$$ means the *N*-by-1
vector with (*i*,1)-component

$$\sum _{j=1}^{N}\left(\mathrm{cos}(\alpha ){c}_{i,j,1,1}\frac{\partial}{\partial x}+\mathrm{cos}(\alpha ){c}_{i,j,1,2}\frac{\partial}{\partial y}+\mathrm{sin}(\alpha ){c}_{i,j,2,1}\frac{\partial}{\partial x}+\mathrm{sin}(\alpha ){c}_{i,j,2,2}\frac{\partial}{\partial y}\right)\text{\hspace{0.17em}}}{u}_{j$$

where the outward normal vector of the boundary $$n=\left(\mathrm{cos}(\alpha ),\mathrm{sin}(\alpha )\right)$$.

For 3-D systems, the notation $$n\text{\hspace{0.17em}}\xb7\text{\hspace{0.17em}}\left(c\otimes \nabla u\right)$$ means
the *N*-by-1 vector with (*i*,1)-component

$$\begin{array}{l}{\displaystyle \sum _{j=1}^{N}\left(\mathrm{sin}\left(\phi \right)\mathrm{cos}\left(\theta \right){c}_{i,j,1,1}\frac{\partial}{\partial x}+\mathrm{sin}\left(\phi \right)\mathrm{cos}\left(\theta \right){c}_{i,j,1,2}\frac{\partial}{\partial y}+\mathrm{sin}\left(\phi \right)\mathrm{cos}\left(\theta \right){c}_{i,j,1,3}\frac{\partial}{\partial z}\right){u}_{j}}\\ +{\displaystyle \sum _{j=1}^{N}\left(\mathrm{sin}\left(\phi \right)\mathrm{sin}\left(\theta \right){c}_{i,j,2,1}\frac{\partial}{\partial x}+\mathrm{sin}\left(\phi \right)\mathrm{sin}\left(\theta \right){c}_{i,j,2,2}\frac{\partial}{\partial y}+\mathrm{sin}\left(\phi \right)\mathrm{sin}\left(\theta \right){c}_{i,j,2,3}\frac{\partial}{\partial z}\right){u}_{j}}\\ +{\displaystyle \sum _{j=1}^{N}\left(\mathrm{cos}\left(\theta \right){c}_{i,j,3,1}\frac{\partial}{\partial x}+\mathrm{cos}\left(\theta \right){c}_{i,j,3,2}\frac{\partial}{\partial y}+\mathrm{cos}\left(\theta \right){c}_{i,j,3,3}\frac{\partial}{\partial z}\right){u}_{j}}\end{array}$$

where the outward normal vector of the boundary $$n=\left(\mathrm{sin}(\phi )\mathrm{cos}(\theta ),\mathrm{sin}(\phi )\mathrm{sin}(\theta ),\mathrm{cos}(\phi )\right)$$. For each edge or face segment,
there are a total of *N* boundary conditions.

Suppose that you have a PDE model named `model`

,
and edges or faces `[e1,e2,e3]`

where the first component
of the solution *u* must satisfy the Neumann boundary
condition with `q = 2`

and `g = 3`

,
and the second component must satisfy the Neumann boundary condition
with `q = 4`

and `g = 5`

. Specify
this boundary condition as follows.

Q = [2 0; 0 4]; G = [3;5]; % For 3-D geometry: applyBoundaryCondition(model,'neumann','Face',[e1,e2,e3],'q',Q,'g',G); % For 2-D geometry: applyBoundaryCondition(model,'neumann','Edge',[e1,e2,e3],'q',Q,'g',G);

The

`'g'`

parameter must be a numeric vector of length*N*. If you do not specify`'g'`

,`applyBoundaryCondition`

sets the values to`0`

.The

`'q'`

parameter can be an*N*-by-*N*numeric matrix or a vector of length*N*^{2}corresponding to the linear indexing form of the*N*-by-*N*matrix. For details about the linear indexing form, see Array Indexing (MATLAB). If you do not specify`'q'`

,`applyBoundaryCondition`

sets the values to`0`

.

If some equations in your system of PDEs must satisfy the Dirichlet
boundary condition and some must satisfy the Neumann boundary condition
for the same geometric region, use the `'mixed'`

parameter
to apply boundary conditions in one call. Note that `applyBoundaryCondition`

uses
the default Neumann boundary condition with `g = 0`

and ```
q
= 0
```

for equations for which you do not explicitly specify
a boundary condition.

Suppose that you have a PDE model named `model`

,
and edge or face labels `[e1,e2,e3]`

where the first
component of the solution *u* must equal `11`

,
the second component must equal `22`

, and the third
component must satisfy the Neumann boundary condition with ```
q
= 3
```

and `g = 4`

. Express this boundary
condition as follows.

Q = [0 0 0; 0 0 0; 0 0 3]; G = [0;0;4]; % For 3-D geometry: applyBoundaryCondition(model,'mixed','Face',[e1,e2,e3],... 'u',[11,22],'EquationIndex',[1,2],... 'q',Q,'g',G); % For 2-D geometry: applyBoundaryCondition(model,'mixed',... 'Edge',[e1,e2,e3],'u',[11,22],... 'EquationIndex',[1,2],'q',Q,'g',G);

Suppose that you have a PDE model named `model`

,
and edges or faces `[e1,e2,e3]`

where the first component
of the solution *u* must satisfy the Dirichlet boundary
condition 2*u _{1}* = 3, the second
component must satisfy the Neumann boundary condition with

```
q
= 4
```

and `g = 5`

, and the third component
must satisfy the Neumann boundary condition with `q = 6`

and ```
g
= 7
```

. Express this boundary condition as follows.h = [2 0 0; 0 0 0; 0 0 0]; r = [3;0;0]; Q = [0 0 0; 0 4 0; 0 0 6]; G = [0;5;7]; % For 3-D geometry: applyBoundaryCondition(model,'mixed', ... 'Face',[e1,e2,e3], ... 'h',h,'r',r,'q',Q,'g',G); % For 2-D geometry: applyBoundaryCondition(model,'mixed', ... 'Edge',[e1,e2,e3], ... 'h',h,'r',r,'q',Q,'g',G);

Use functions to express nonconstant boundary conditions.

applyBoundaryCondition(model,'dirichlet', ... 'Edge',1, ... 'r',@myrfun); applyBoundaryCondition(model,'neumann', ... 'Face',2, ... 'g',@mygfun,'q',@myqfun); applyBoundaryCondition(model,'mixed', ... 'Edge',[3,4], ... 'u',@myufun, ... 'EquationIndex',[2,3]);

Each function must have the following syntax.

`function bcMatrix = myfun(location,state)`

Partial Differential Equation Toolbox™ solvers pass the `location`

and
`state`

data to your function.

`location`

— A structure containing the following fields. If you pass a name-value pair to`applyBoundaryCondition`

with`Vectorized`

set to`'on'`

, then`location`

can contain several evaluation points. If you do not set`Vectorized`

or use`Vectorized,'off'`

, then solvers pass just one evaluation point in each call.`location.x`

— The*x*-coordinate of the point or points`location.y`

— The*y*-coordinate of the point or points`location.z`

— For 3-D geometry, the*z*-coordinate of the point or points

Furthermore, if there are Neumann conditions, then solvers pass the following data in the

`location`

structure.`location.nx`

—*x*-component of the normal vector at the evaluation point or points`location.ny`

—*y*-component of the normal vector at the evaluation point or points`location.nz`

— For 3-D geometry,*z*-component of the normal vector at the evaluation point or points

`state`

— For transient or nonlinear problems.`state.u`

contains the solution vector at evaluation points.`state.u`

is an*N*-by-*M*matrix, where each column corresponds to one evaluation point, and*M*is the number of evaluation points.`state.time`

contains the time at evaluation points.`state.time`

is a scalar.

Your function returns `bcMatrix`

. This matrix
has the following form, depending on the boundary condition type.

`'u'`

—*N1*-by-*M*matrix, where each column corresponds to one evaluation point, and*M*is the number of evaluation points.*N1*is the length of the`'EquationIndex'`

argument. If there is no`'EquationIndex'`

argument, then*N1*=*N*.`'r'`

or`'g'`

—*N*-by-*M*matrix, where each column corresponds to one evaluation point, and*M*is the number of evaluation points.`'h'`

or`'q'`

—*N*^{2}-by-*M*matrix, where each column corresponds to one evaluation point via linear indexing of the underlying*N*-by-*N*matrix, and*M*is the number of evaluation points. Alternatively, an*N*-by-*N*-by-*M*array, where each evaluation point is an*N*-by-*N*matrix. For details about linear indexing, see Array Indexing (MATLAB).

If boundary conditions depend on `state.u`

or `state.time`

,
ensure that your function returns a matrix of `NaN`

of
the correct size when `state.u`

or `state.time`

are `NaN`

.
Solvers check whether a problem is nonlinear or time-dependent by
passing `NaN`

state values, and looking for returned `NaN`

values.