NAG Library Manual, Mark 29.2
```!   D03EEF Example Program Text
!   Mark 29.2 Release. NAG Copyright 2023.

Module d03eefe_mod

!     D03EEF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Accessibility Statements ..
Private
Public                           :: bndy, fexact, pdef
!     .. Parameters ..
Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
Integer, Parameter, Public       :: nin = 5, nout = 6
Contains
Subroutine pdef(x,y,alpha,beta,gamma,delta,epslon,phi,psi)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Out) :: alpha, beta, delta, epslon, gamma, &
phi, psi
Real (Kind=nag_wp), Intent (In) :: x, y
!       .. Intrinsic Procedures ..
Intrinsic                      :: cos, sin
!       .. Executable Statements ..
alpha = one
beta = zero
gamma = one
delta = 50.0_nag_wp
epslon = 50.0_nag_wp
phi = zero

psi = sin(x)*((-alpha-gamma+phi)*sin(y)+epslon*cos(y)) +               &
cos(x)*(delta*sin(y)+beta*cos(y))

Return
End Subroutine pdef
Subroutine bndy(x,y,a,b,c,ibnd)

!       .. Parameters ..
Integer, Parameter             :: bottom = 0, left = 3, right = 1,     &
top = 2
!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Out) :: a, b, c
Real (Kind=nag_wp), Intent (In) :: x, y
Integer, Intent (In)           :: ibnd
!       .. Intrinsic Procedures ..
Intrinsic                      :: sin
!       .. Executable Statements ..
If (ibnd==top .Or. ibnd==right) Then

!         Solution prescribed

a = one
b = zero
c = sin(x)*sin(y)
Else If (ibnd==bottom) Then

!         Derivative prescribed

a = zero
b = one
c = -sin(x)
Else If (ibnd==left) Then

!         Derivative prescribed

a = zero
b = one
c = -sin(y)
End If

Return
End Subroutine bndy
Function fexact(x,y)

!       .. Function Return Value ..
Real (Kind=nag_wp)             :: fexact
!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: x, y
!       .. Intrinsic Procedures ..
Intrinsic                      :: sin
!       .. Executable Statements ..
fexact = sin(x)*sin(y)
Return
End Function fexact
End Module d03eefe_mod

Program d03eefe

!     D03EEF Example Main Program

!     .. Use Statements ..
Use d03eefe_mod, Only: bndy, fexact, nin, nout, pdef, zero
Use nag_library, Only: d03edf, d03eef, nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: acc, hx, hy, rmserr, xmax, xmin, xx, &
ymax, ymin, yy
Integer                          :: i, icase, ifail, iout, ix, j, lda,   &
levels, maxit, ngx, ngxy, ngy, numit
Character (7)                    :: scheme
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: a(:,:), rhs(:), u(:), ub(:), us(:),  &
x(:), y(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: real, sqrt
!     .. Executable Statements ..
Write (nout,*) 'D03EEF Example Program Results'
Write (nout,*)
Flush (nout)

!     Skip heading in data file
ngx = 2**levels + 1
ngy = ngx
lda = 4*(ngx+1)*(ngy+1)/3
ngxy = ngx*ngy
Allocate (a(lda,7),rhs(lda),u(lda),ub(ngxy),us(lda),x(ngxy),y(ngxy))
hx = (xmax-xmin)/real(ngx-1,kind=nag_wp)
Do i = 1, ngx
xx = xmin + real(i-1,kind=nag_wp)*hx
x(i:ngxy:ngx) = xx
End Do
hy = (ymax-ymin)/real(ngy-1,kind=nag_wp)
Do j = 1, ngy
yy = ymin + real(j-1,kind=nag_wp)*hy
y((j-1)*ngx+1:j*ngx) = yy
End Do
!     ** set iout > 2 to obtain intermediate output from D03EDF **
iout = 0

cases: Do icase = 1, 2

Select Case (icase)
Case (1)
!         Central differences
scheme = 'Central'
Case (2)
!         Upwind differences
scheme = 'Upwind'
End Select
!       Discretize the equations
!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = -1
Call d03eef(xmin,xmax,ymin,ymax,pdef,bndy,ngx,ngy,lda,a,rhs,scheme,    &
ifail)

If (ifail<0) Then
Write (nout,99995) ifail
Exit cases
End If

!       Set the initial guess to zero
ub(1:ngxy) = zero

!       Solve the equations
ifail = 0
Call d03edf(ngx,ngy,lda,a,rhs,ub,maxit,acc,us,u,iout,numit,ifail)

!       Print out the solution
Write (nout,*)
Write (nout,*) 'Exact solution above computed solution'
Write (nout,*)
Write (nout,99998) '  I/J', (i,i=1,ngx)
rmserr = zero
Do j = ngy, 1, -1
ix = (j-1)*ngx
Write (nout,*)
Write (nout,99999) j, (fexact(x(ix+i),y(ix+i)),i=1,ngx)
Write (nout,99999) j, u(ix+1:ix+ngx)
Do i = 1, ngx
rmserr = rmserr + (fexact(x(ix+i),y(ix+i))-u(ix+i))**2
End Do
End Do
rmserr = sqrt(rmserr/real(ngxy,kind=nag_wp))
Write (nout,*)
Write (nout,99997) 'Number of Iterations = ', numit
Write (nout,99996) 'RMS Error = ', rmserr
End Do cases

99999 Format (1X,I3,2X,10F7.3,:,/,(6X,10F7.3))
99998 Format (1X,A,10I7,:,/,(6X,10I7))
99997 Format (1X,A,I3)
99996 Format (1X,A,1P,E10.2)
99995 Format (1X,' ** D03EEF returned with IFAIL = ',I5)
End Program d03eefe
```