```    Program d03fafe

!     D03FAF Example Program Text

!     Mark 27.0 Release. NAG Copyright 2019.

!     .. Use Statements ..
Use nag_library, Only: d03faf, nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: dx, dy, dz, error, lambda, pertrb,   &
t, xf, xs, yf, ys, yy, zf, zs, zz
Integer                          :: i, ifail, j, k, l, lbdcnd, ldf,      &
ldf2, lwrk, m, mbdcnd, n, nbdcnd
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: bdxf(:,:), bdxs(:,:), bdyf(:,:),     &
bdys(:,:), bdzf(:,:), bdzs(:,:),     &
f(:,:,:), x(:), y(:), z(:)
Real (Kind=nag_wp)               :: w(0)
!     .. Intrinsic Procedures ..
Intrinsic                        :: abs, cos, real, sin
!     .. Executable Statements ..
Write (nout,*) 'D03FAF Example Program Results'
!     Skip heading in data file
ldf = l + 1
ldf2 = m + 1
lwrk = 0
Allocate (bdxf(ldf2,n+1),bdxs(ldf2,n+1),bdyf(ldf,n+1),bdys(ldf,n+1),     &
bdzf(ldf,m+1),bdzs(ldf,m+1),f(ldf,ldf2,n+1),x(l+1),y(m+1),z(n+1))

!     Define the grid points for later use.

dx = (xf-xs)/real(l,kind=nag_wp)
Do i = 1, l + 1
x(i) = xs + real(i-1,kind=nag_wp)*dx
End Do
dy = (yf-ys)/real(m,kind=nag_wp)
Do j = 1, m + 1
y(j) = ys + real(j-1,kind=nag_wp)*dy
End Do
dz = (zf-zs)/real(n,kind=nag_wp)
Do k = 1, n + 1
z(k) = zs + real(k-1,kind=nag_wp)*dz
End Do

!     Define the array of derivative boundary values.

Do j = 1, m + 1
yy = sin(y(j))
bdzf(1:l+1,j) = -yy*x(1:l+1)**4
End Do

!     Note that for this example all other boundary arrays are
!     dummy variables.

!     We define the function boundary values in the F array.

f(1,1:m+1,1:n+1) = 0.0_nag_wp
Do k = 1, n + 1
zz = cos(z(k))
Do j = 1, m + 1
f(l+1,j,k) = zz*sin(y(j))
End Do
End Do
f(1:l+1,1:m+1,1) = -bdzf(1:l+1,1:m+1)

!     Define the values of the right hand side of the Helmholtz
!     equation.

Do k = 2, n + 1
zz = 4.0_nag_wp*cos(z(k))
Do j = 1, m + 1
yy = sin(y(j))*zz
Do i = 2, l
f(i,j,k) = yy*x(i)**2*(3.0_nag_wp-x(i)**2)
End Do
End Do
End Do

!     Call D03FAF to generate and solve the finite difference equation.
!     ifail: behaviour on error exit
!            =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d03faf(xs,xf,l,lbdcnd,bdxs,bdxf,ys,yf,m,mbdcnd,bdys,bdyf,zs,zf,n,   &
nbdcnd,bdzs,bdzf,lambda,ldf,ldf2,f,pertrb,w,lwrk,ifail)

!     Compute discretization error.  The exact solution to the
!     problem is

!     U(X,Y,Z) = X**4*SIN(Y)*COS(Z)

error = 0.0_nag_wp
Do k = 1, n + 1
zz = cos(z(k))
Do j = 1, m + 1
yy = sin(y(j))*zz
Do i = 1, l + 1
t = abs(f(i,j,k)-yy*x(i)**4)
If (t>error) Then
error = t
End If
End Do
End Do
End Do
Write (nout,*)
Write (nout,99999) error

99999 Format (1X,'Maximum component of discretization error =',1P,E13.6)
End Program d03fafe
```