Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_pde_3d_ellip_helmholtz (d03fa)

## Purpose

nag_pde_3d_ellip_helmholtz (d03fa) solves the Helmholtz equation in Cartesian coordinates in three dimensions using the standard seven-point finite difference approximation. This function is designed to be particularly efficient on vector processors.

## Syntax

[f, pertrb, ifail] = d03fa(xs, xf, l, lbdcnd, bdxs, bdxf, ys, yf, m, mbdcnd, bdys, bdyf, zs, zf, n, nbdcnd, bdzs, bdzf, lambda, f, 'lwrk', lwrk)
[f, pertrb, ifail] = nag_pde_3d_ellip_helmholtz(xs, xf, l, lbdcnd, bdxs, bdxf, ys, yf, m, mbdcnd, bdys, bdyf, zs, zf, n, nbdcnd, bdzs, bdzf, lambda, f, 'lwrk', lwrk)
Note: the interface to this routine has changed since earlier releases of the toolbox:
 At Mark 23: lwrk was added as an optional parameter

## Description

nag_pde_3d_ellip_helmholtz (d03fa) solves the three-dimensional Helmholtz equation in Cartesian coordinates:
 $∂2u ∂x2 + ∂2u ∂y2 + ∂2u ∂z2 +λu=fx,y,z.$
This function forms the system of linear equations resulting from the standard seven-point finite difference equations, and then solves the system using a method based on the fast Fourier transform (FFT) described by Swarztrauber (1984). This function is based on the function HW3CRT from FISHPACK (see Swarztrauber and Sweet (1979)).
More precisely, the function replaces all the second derivatives by second-order central difference approximations, resulting in a block tridiagonal system of linear equations. The equations are modified to allow for the prescribed boundary conditions. Either the solution or the derivative of the solution may be specified on any of the boundaries, or the solution may be specified to be periodic in any of the three dimensions. By taking the discrete Fourier transform in the $x$- and $y$-directions, the equations are reduced to sets of tridiagonal systems of equations. The Fourier transforms required are computed using the multiple FFT functions found in Chapter C06.

## References

Swarztrauber P N (1984) Fast Poisson solvers Studies in Numerical Analysis (ed G H Golub) Mathematical Association of America
Swarztrauber P N and Sweet R A (1979) Efficient Fortran subprograms for the solution of separable elliptic partial differential equations ACM Trans. Math. Software 5 352–364

## Parameters

### Compulsory Input Parameters

1:     $\mathrm{xs}$ – double scalar
The lower bound of the range of $x$, i.e., ${\mathbf{xs}}\le x\le {\mathbf{xf}}$.
Constraint: ${\mathbf{xs}}<{\mathbf{xf}}$.
2:     $\mathrm{xf}$ – double scalar
The upper bound of the range of $x$, i.e., ${\mathbf{xs}}\le x\le {\mathbf{xf}}$.
Constraint: ${\mathbf{xs}}<{\mathbf{xf}}$.
3:     $\mathrm{l}$int64int32nag_int scalar
The number of panels into which the interval (xs,xf) is subdivided. Hence, there will be ${\mathbf{l}}+1$ grid points in the $x$-direction given by ${x}_{\mathit{i}}={\mathbf{xs}}+\left(\mathit{i}-1\right)×\delta x$, for $\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$, where $\delta x=\left({\mathbf{xf}}-{\mathbf{xs}}\right)/{\mathbf{l}}$ is the panel width.
Constraint: ${\mathbf{l}}\ge 5$.
4:     $\mathrm{lbdcnd}$int64int32nag_int scalar
Indicates the type of boundary conditions at $x={\mathbf{xs}}$ and $x={\mathbf{xf}}$.
${\mathbf{lbdcnd}}=0$
If the solution is periodic in $x$, i.e., $u\left({\mathbf{xs}},y,z\right)=u\left({\mathbf{xf}},y,z\right)$.
${\mathbf{lbdcnd}}=1$
If the solution is specified at $x={\mathbf{xs}}$ and $x={\mathbf{xf}}$.
${\mathbf{lbdcnd}}=2$
If the solution is specified at $x={\mathbf{xs}}$ and the derivative of the solution with respect to $x$ is specified at $x={\mathbf{xf}}$.
${\mathbf{lbdcnd}}=3$
If the derivative of the solution with respect to $x$ is specified at $x={\mathbf{xs}}$ and $x={\mathbf{xf}}$.
${\mathbf{lbdcnd}}=4$
If the derivative of the solution with respect to $x$ is specified at $x={\mathbf{xs}}$ and the solution is specified at $x={\mathbf{xf}}$.
Constraint: $0\le {\mathbf{lbdcnd}}\le 4$.
5:     $\mathrm{bdxs}\left(\mathit{ldf2},{\mathbf{n}}+1\right)$ – double array
ldf2, the first dimension of the array, must satisfy the constraint $\mathit{ldf2}\ge {\mathbf{m}}+1$.
The values of the derivative of the solution with respect to $x$ at $x={\mathbf{xs}}$. When ${\mathbf{lbdcnd}}=3$ or $4$, ${\mathbf{bdxs}}\left(\mathit{j},\mathit{k}\right)={u}_{x}\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$, for $\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$ and $\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$.
When lbdcnd has any other value, bdxs is not referenced.
6:     $\mathrm{bdxf}\left(\mathit{ldf2},{\mathbf{n}}+1\right)$ – double array
ldf2, the first dimension of the array, must satisfy the constraint $\mathit{ldf2}\ge {\mathbf{m}}+1$.
The values of the derivative of the solution with respect to $x$ at $x={\mathbf{xf}}$. When ${\mathbf{lbdcnd}}=2$ or $3$, ${\mathbf{bdxf}}\left(\mathit{j},\mathit{k}\right)={u}_{x}\left({\mathbf{xf}},{y}_{\mathit{j}},{z}_{\mathit{k}}\right)$, for $\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$ and $\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$.
When lbdcnd has any other value, bdxf is not referenced.
7:     $\mathrm{ys}$ – double scalar
The lower bound of the range of $y$, i.e., ${\mathbf{ys}}\le y\le {\mathbf{yf}}$.
Constraint: ${\mathbf{ys}}<{\mathbf{yf}}$.
8:     $\mathrm{yf}$ – double scalar
The upper bound of the range of $y$, i.e., ${\mathbf{ys}}\le y\le {\mathbf{yf}}$.
Constraint: ${\mathbf{ys}}<{\mathbf{yf}}$.
9:     $\mathrm{m}$int64int32nag_int scalar
The number of panels into which the interval (ys,yf) is subdivided. Hence, there will be ${\mathbf{m}}+1$ grid points in the $y$-direction given by ${y}_{\mathit{j}}={\mathbf{ys}}+\left(\mathit{j}-1\right)×\delta y$, for $\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$, where $\delta y=\left({\mathbf{yf}}-{\mathbf{ys}}\right)/{\mathbf{m}}$ is the panel width.
Constraint: ${\mathbf{m}}\ge 5$.
10:   $\mathrm{mbdcnd}$int64int32nag_int scalar
Indicates the type of boundary conditions at $y={\mathbf{ys}}$ and $y={\mathbf{yf}}$.
${\mathbf{mbdcnd}}=0$
If the solution is periodic in $y$, i.e., $u\left(x,{\mathbf{yf}},z\right)=u\left(x,{\mathbf{ys}},z\right)$.
${\mathbf{mbdcnd}}=1$
If the solution is specified at $y={\mathbf{ys}}$ and $y={\mathbf{yf}}$.
${\mathbf{mbdcnd}}=2$
If the solution is specified at $y={\mathbf{ys}}$ and the derivative of the solution with respect to $y$ is specified at $y={\mathbf{yf}}$.
${\mathbf{mbdcnd}}=3$
If the derivative of the solution with respect to $y$ is specified at $y={\mathbf{ys}}$ and $y={\mathbf{yf}}$.
${\mathbf{mbdcnd}}=4$
If the derivative of the solution with respect to $y$ is specified at $y={\mathbf{ys}}$ and the solution is specified at $y={\mathbf{yf}}$.
Constraint: $0\le {\mathbf{mbdcnd}}\le 4$.
11:   $\mathrm{bdys}\left(\mathit{ldf},{\mathbf{n}}+1\right)$ – double array
ldf, the first dimension of the array, must satisfy the constraint $\mathit{ldf}\ge {\mathbf{l}}+1$.
The values of the derivative of the solution with respect to $y$ at $y={\mathbf{ys}}$. When ${\mathbf{mbdcnd}}=3$ or $4$, ${\mathbf{bdys}}\left(\mathit{i},\mathit{k}\right)={u}_{y}\left({x}_{\mathit{i}},{\mathbf{ys}},{z}_{\mathit{k}}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$ and $\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$.
When mbdcnd has any other value, bdys is not referenced.
12:   $\mathrm{bdyf}\left(\mathit{ldf},{\mathbf{n}}+1\right)$ – double array
ldf, the first dimension of the array, must satisfy the constraint $\mathit{ldf}\ge {\mathbf{l}}+1$.
The values of the derivative of the solution with respect to $y$ at $y={\mathbf{yf}}$. When ${\mathbf{mbdcnd}}=2$ or $3$, ${\mathbf{bdyf}}\left(\mathit{i},\mathit{k}\right)={u}_{y}\left({x}_{\mathit{i}},{\mathbf{yf}},{z}_{\mathit{k}}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$ and $\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$.
When mbdcnd has any other value, bdyf is not referenced.
13:   $\mathrm{zs}$ – double scalar
The lower bound of the range of $z$, i.e., ${\mathbf{zs}}\le z\le {\mathbf{zf}}$.
Constraint: ${\mathbf{zs}}<{\mathbf{zf}}$.
14:   $\mathrm{zf}$ – double scalar
The upper bound of the range of $z$, i.e., ${\mathbf{zs}}\le z\le {\mathbf{zf}}$.
Constraint: ${\mathbf{zs}}<{\mathbf{zf}}$.
15:   $\mathrm{n}$int64int32nag_int scalar
The number of panels into which the interval (zs,zf) is subdivided. Hence, there will be ${\mathbf{n}}+1$ grid points in the $z$-direction given by ${z}_{\mathit{k}}={\mathbf{zs}}+\left(\mathit{k}-1\right)×\delta z$, for $\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$, where $\delta z=\left({\mathbf{zf}}-{\mathbf{zs}}\right)/{\mathbf{n}}$ is the panel width.
Constraint: ${\mathbf{n}}\ge 5$.
16:   $\mathrm{nbdcnd}$int64int32nag_int scalar
Specifies the type of boundary conditions at $z={\mathbf{zs}}$ and $z={\mathbf{zf}}$.
${\mathbf{nbdcnd}}=0$
if the solution is periodic in $z$, i.e., $u\left(x,y,{\mathbf{zf}}\right)=u\left(x,y,{\mathbf{zs}}\right)$.
${\mathbf{nbdcnd}}=1$
if the solution is specified at $z={\mathbf{zs}}$ and $z={\mathbf{zf}}$.
${\mathbf{nbdcnd}}=2$
if the solution is specified at $z={\mathbf{zs}}$ and the derivative of the solution with respect to $z$ is specified at $z={\mathbf{zf}}$.
${\mathbf{nbdcnd}}=3$
if the derivative of the solution with respect to $z$ is specified at $z={\mathbf{zs}}$ and $z={\mathbf{zf}}$.
${\mathbf{nbdcnd}}=4$
if the derivative of the solution with respect to $z$ is specified at $z={\mathbf{zs}}$ and the solution is specified at $z={\mathbf{zf}}$.
Constraint: $0\le {\mathbf{nbdcnd}}\le 4$.
17:   $\mathrm{bdzs}\left(\mathit{ldf},{\mathbf{m}}+1\right)$ – double array
ldf, the first dimension of the array, must satisfy the constraint $\mathit{ldf}\ge {\mathbf{l}}+1$.
The values of the derivative of the solution with respect to $z$ at $z={\mathbf{zs}}$. When ${\mathbf{nbdcnd}}=3$ or $4$, ${\mathbf{bdzs}}\left(\mathit{i},\mathit{j}\right)={u}_{z}\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$ and $\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$.
When nbdcnd has any other value, bdzs is not referenced.
18:   $\mathrm{bdzf}\left(\mathit{ldf},{\mathbf{m}}+1\right)$ – double array
ldf, the first dimension of the array, must satisfy the constraint $\mathit{ldf}\ge {\mathbf{l}}+1$.
The values of the derivative of the solution with respect to $z$ at $z={\mathbf{zf}}$. When ${\mathbf{nbdcnd}}=2$ or $3$, ${\mathbf{bdzf}}\left(\mathit{i},\mathit{j}\right)={u}_{z}\left({x}_{i},{y}_{j},{\mathbf{zf}}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$ and $\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$.
When nbdcnd has any other value, bdzf is not referenced.
19:   $\mathrm{lambda}$ – double scalar
The constant $\lambda$ in the Helmholtz equation. For certain positive values of $\lambda$ a solution to the differential equation may not exist, and close to these values the solution of the discretized problem will be extremely ill-conditioned. If $\lambda >0$, then nag_pde_3d_ellip_helmholtz (d03fa) will set ${\mathbf{ifail}}={\mathbf{3}}$, but will still attempt to find a solution. However, since in general the values of $\lambda$ for which no solution exists cannot be predicted a priori, you are advised to treat any results computed with $\lambda >0$ with great caution.
20:   $\mathrm{f}\left(\mathit{ldf},\mathit{ldf2},{\mathbf{n}}+1\right)$ – double array
ldf, the first dimension of the array, must satisfy the constraint $\mathit{ldf}\ge {\mathbf{l}}+1$.
The values of the right-side of the Helmholtz equation and boundary values (if any).
 $fijk = fxi,yj,zk i=2,3,…,l, j=2,3,…,m ​ and ​ k=2,3,…,n.$
On the boundaries f is defined by
 lbdcnd ${\mathbf{f}}\left(1,j,k\right)$ ${\mathbf{f}}\left({\mathbf{l}}+1,j,k\right)$ 0 $f\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ $f\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ 1 $u\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ $u\left({\mathbf{xf}},{y}_{j},{z}_{k}\right)$ 2 $u\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ $f\left({\mathbf{xf}},{y}_{j},{z}_{k}\right)$ $j=1,2,\dots ,{\mathbf{m}}+1$ 3 $f\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ $f\left({\mathbf{xf}},{y}_{j},{z}_{k}\right)$ $k=1,2,\dots ,{\mathbf{n}}+1$ 4 $f\left({\mathbf{xs}},{y}_{j},{z}_{k}\right)$ $u\left({\mathbf{xf}},{y}_{j},{z}_{k}\right)$ mbdcnd ${\mathbf{f}}\left(i,1,k\right)$ ${\mathbf{f}}\left(i,{\mathbf{m}}+1,k\right)$ 0 $f\left({x}_{i},{\mathbf{ys}},{z}_{k}\right)$ $f\left({x}_{i},{\mathbf{ys}},{z}_{k}\right)$ 1 $u\left({\mathbf{ys}},{x}_{i},{z}_{k}\right)$ $u\left({\mathbf{yf}},{x}_{i},{z}_{k}\right)$ 2 $u\left({x}_{i},{\mathbf{ys}},{z}_{k}\right)$ $f\left({x}_{i},{\mathbf{yf}},{z}_{k}\right)$ $i=1,2,\dots ,{\mathbf{l}}+1$ 3 $f\left({x}_{i},{\mathbf{ys}},{z}_{k}\right)$ $f\left({x}_{i},{\mathbf{yf}},{z}_{k}\right)$ $k=1,2,\dots ,{\mathbf{n}}+1$ 4 $f\left({x}_{i},{\mathbf{ys}},{z}_{k}\right)$ $u\left({x}_{i},{\mathbf{yf}},{z}_{k}\right)$ nbdcnd ${\mathbf{f}}\left(i,j,1\right)$ ${\mathbf{f}}\left(i,j,{\mathbf{n}}+1\right)$ 0 $f\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ $f\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ 1 $u\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ $u\left({x}_{i},{y}_{j},{\mathbf{zf}}\right)$ 2 $u\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ $f\left({x}_{i},{y}_{j},{\mathbf{zf}}\right)$ $i=1,2,\dots ,{\mathbf{l}}+1$ 3 $f\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ $f\left({x}_{i},{y}_{j},{\mathbf{zf}}\right)$ $j=1,2,\dots ,{\mathbf{m}}+1$ 4 $f\left({x}_{i},{y}_{j},{\mathbf{zs}}\right)$ $u\left({x}_{i},{y}_{j},{\mathbf{zf}}\right)$
Note: if the table calls for both the solution $u$ and the right-hand side $f$ on a boundary, then the solution must be specified.

### Optional Input Parameters

1:     $\mathrm{lwrk}$int64int32nag_int scalar
Default: $2×\left({\mathbf{n}}+1\right)×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{l}},{\mathbf{m}}\right)+3×{\mathbf{l}}+3×{\mathbf{m}}+4×{\mathbf{n}}+6$
The dimension of the array w. $2×\left({\mathbf{n}}+1\right)×\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{l}},{\mathbf{m}}\right)+3×{\mathbf{l}}+3×{\mathbf{m}}+4×{\mathbf{n}}+6$ is an upper bound on the required size of w. If lwrk is too small, the function exits with ${\mathbf{ifail}}={\mathbf{2}}$, and if on entry ${\mathbf{ifail}}=0$ or $-1$, a message is output giving the exact value of lwrk required to solve the current problem.

### Output Parameters

1:     $\mathrm{f}\left(\mathit{ldf},\mathit{ldf2},{\mathbf{n}}+1\right)$ – double array
$\mathit{ldf2}={\mathbf{m}}+1$.
Contains the solution $u\left(\mathit{i},\mathit{j},\mathit{k}\right)$ of the finite difference approximation for the grid point $\left({x}_{\mathit{i}},{y}_{\mathit{j}},{z}_{\mathit{k}}\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{l}}+1$, $\mathit{j}=1,2,\dots ,{\mathbf{m}}+1$ and $\mathit{k}=1,2,\dots ,{\mathbf{n}}+1$.
2:     $\mathrm{pertrb}$ – double scalar
${\mathbf{pertrb}}=0$, unless a solution to Poisson's equation $\left(\lambda =0\right)$ is required with a combination of periodic or derivative boundary conditions (lbdcnd, mbdcnd and ${\mathbf{nbdcnd}}=0$ or $3$). In this case a solution may not exist. pertrb is a constant, calculated and subtracted from the array f, which ensures that a solution exists. nag_pde_3d_ellip_helmholtz (d03fa) then computes this solution, which is a least squares solution to the original approximation. This solution is not unique and is unnormalized. The value of pertrb should be small compared to the right-hand side f, otherwise a solution has been obtained to an essentially different problem. This comparison should always be made to ensure that a meaningful solution has been obtained.
3:     $\mathrm{ifail}$int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see Error Indicators and Warnings).

## Error Indicators and Warnings

Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

${\mathbf{ifail}}=1$
 On entry, ${\mathbf{xs}}\ge {\mathbf{xf}}$, or ${\mathbf{l}}<5$, or ${\mathbf{lbdcnd}}<0$, or ${\mathbf{lbdcnd}}>4$, or ${\mathbf{ys}}\ge {\mathbf{yf}}$, or ${\mathbf{m}}<5$, or ${\mathbf{mbdcnd}}<0$, or ${\mathbf{mbdcnd}}>4$, or ${\mathbf{zs}}\ge {\mathbf{zf}}$, or ${\mathbf{n}}<5$, or ${\mathbf{nbdcnd}}<0$, or ${\mathbf{nbdcnd}}>4$, or $\mathit{ldf}<{\mathbf{l}}+1$, or $\mathit{ldf2}<{\mathbf{m}}+1$.
${\mathbf{ifail}}=2$
 On entry, lwrk is too small.
W  ${\mathbf{ifail}}=3$
 On entry, $\lambda >0$.
${\mathbf{ifail}}=-99$
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

Not applicable.

The execution time is roughly proportional to ${\mathbf{l}}×{\mathbf{m}}×{\mathbf{n}}×\left({\mathrm{log}}_{2}{\mathbf{l}}+{\mathrm{log}}_{2}{\mathbf{m}}+5\right)$, but also depends on input arguments lbdcnd and mbdcnd.

## Example

This example solves the Helmholz equation
 $∂2u ∂x2 + ∂2u ∂y2 + ∂2u ∂z2 +λ u=fx,y,z$
for $\left(x,y,z\right)\in \left[0,1\right]×\left[0,2\pi \right]×\left[0,\frac{\pi }{2}\right]$, where $\lambda =-2$, and $f\left(x,y,z\right)$ is derived from the exact solution
 $ux,y,z=x4sin⁡ycos⁡z.$
The equation is subject to the following boundary conditions, again derived from the exact solution given above.
• $u\left(0,y,z\right)$ and $u\left(1,y,z\right)$ are prescribed (i.e., ${\mathbf{lbdcnd}}=1$).
• $u\left(x,0,z\right)=u\left(x,2\pi ,z\right)$ (i.e., ${\mathbf{mbdcnd}}=0$).
• $u\left(x,y,0\right)$ and ${u}_{x}\left(x,y,\frac{\pi }{2}\right)$ are prescribed (i.e., ${\mathbf{nbdcnd}}=2$).
```function d03fa_example

fprintf('d03fa example results\n\n');

xs = 0; xf = 1;    l = 16;
ys = 0; yf = 2*pi; m = 32;
zs = 0; zf = pi/2; n = 20;

bdxs = zeros(m+1, n+1); bdxf = zeros(m+1, n+1);
bdys = zeros(l+1, n+1); bdyf = zeros(l+1, n+1);
bdzs = zeros(l+1, m+1); bdzf = zeros(l+1, m+1);

lbdcnd = int64(1);
mbdcnd = int64(0);
nbdcnd = int64(2);
lambda = -2;

% Define the grid points for later use.
dx = (xf-xs)/l;  x  = [xs:dx:xf];
dy = (yf-ys)/m;  y  = [ys:dy:yf];
dz = (zf-zs)/n;  z  = [zs:dz:zf];

% Define the array of derivative boundary values (z only here).
XY = [x.^4]'*sin(y);
bdzf = -XY;

% Define the f array including boundary conditions.
YZ = [sin(y)]'*cos(z);
for i = 1:l
f(i,:,:) = (4*x(i)^2*(3-x(i)^2))*YZ;
g(i,:,:) = x(i)^4*YZ;
end
g(l+1,:,:) = YZ;
f(1,:,:)   = 0;
f(l+1,:,:) = YZ;
f(:,:,1)   = XY;

[u, pertrb, ifail] = d03fa( ...
xs, xf, int64(l), lbdcnd, bdxs, bdxf, ...
ys, yf, int64(m), mbdcnd, bdys, bdyf, ...
zs, zf, int64(n), nbdcnd, bdzs, bdzf, lambda, f);
% calculate error
maxerr   = max(max(max(abs(u - g))));
fprintf('Maximum error in computed solution = %10.3e\n',maxerr);

fig1 = figure;
[xs,ys,zs] = meshgrid(y,x,z);
slice(xs,ys,zs,u,[1.8 5],[0.95],[0.78]);
axis([y(1) y(end) x(1) x(end) 0 2]);
xlabel('y'); ylabel('x'); zlabel('z');
title('Helholtz Equation solution in a box');
view(111, 26);
h = colorbar;
ylabel(h, 'u(x,y,z)')

```
```d03fa example results

Maximum error in computed solution =  5.177e-04
```