NAG Library Manual, Mark 28.5
```!   D02TGF Example Program Text
!   Mark 28.5 Release. NAG Copyright 2022.

Module d02tgfe_mod

!     D02TGF 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                           :: bdyc, coeff
!     .. Parameters ..
Real (Kind=nag_wp), Parameter, Public :: coeff_tol = 1.0E-5_nag_wp
Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
Real (Kind=nag_wp), Parameter    :: two = 2.0_nag_wp
Integer, Parameter, Public       :: n = 2, nin = 5, nout = 6
!     .. Local Scalars ..
Real (Kind=nag_wp), Public, Save :: x0, x1
Integer, Public, Save            :: k1
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable, Public, Save :: b(:,:)
Integer, Public, Save            :: l(n) = (/1,2/)
Integer, Public, Save            :: m(n) = (/1,2/)
Contains
Subroutine coeff(x,i,a,ia,ia1,rhs)

!       .. Use Statements ..
Use nag_library, Only: e02akf
!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: rhs
Real (Kind=nag_wp), Intent (In) :: x
Integer, Intent (In)           :: i, ia, ia1
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: a(ia,ia1)
!       .. Local Scalars ..
Real (Kind=nag_wp)             :: z1, z2
Integer                        :: ifail
!       .. Executable Statements ..
If (i<=1) Then

!         Evaluate z1, z2 at x using previous coeffs b.
ifail = 0
Call e02akf(k1,x0,x1,b(1,1),1,k1,x,z1,ifail)
Call e02akf(k1,x0,x1,b(1,2),1,k1,x,z2,ifail)

a(1,1) = z2*z2 - one
a(1,2) = two
a(2,1) = two*z1*z2 + one
rhs = two*z1*z2*z2
Else
a(1,2) = -one
a(2,3) = two
End If
Return
End Subroutine coeff
Subroutine bdyc(x,i,j,a,ia,ia1,rhs)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: rhs
Real (Kind=nag_wp), Intent (Out) :: x
Integer, Intent (In)           :: i, ia, ia1, j
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Inout) :: a(ia,ia1)
!       .. Executable Statements ..
x = -one
a(i,j) = one
If (i==2 .And. j==1) Then
rhs = 3.0_nag_wp
End If
Return
End Subroutine bdyc
End Module d02tgfe_mod
Program d02tgfe

!     D02TGF Example Main Program

!     .. Use Statements ..
Use d02tgfe_mod, Only: b, bdyc, coeff, coeff_tol, k1, l, m, n, nin,      &
nout, x0, x1
Use nag_library, Only: d02tgf, e02akf, nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Local Scalars ..
Real (Kind=nag_wp)               :: emax, x, xinc
Integer                          :: i, ia1, ifail, iter, j, k, kp, ldc,  &
liw, lsum, lw, mimax
!     .. Local Arrays ..
Real (Kind=nag_wp), Allocatable  :: c(:,:), w(:), y(:)
Integer, Allocatable             :: iw(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: abs, max, real, sum
!     .. Executable Statements ..
Write (nout,*) 'D02TGF Example Program Results'
!     Skip heading in data file

lsum = sum(l(1:n))
k1 = mimax + 1
ldc = k1
liw = n*k1
lw = 2*(n*kp+lsum)*(n*k1+1) + 7*n*k1
Allocate (b(k1,n),c(ldc,n),w(lw),y(n),iw(liw))

!     Initialize coefficients b(:,:) such that z1 = 0 and z2 = 3.
b(1:k1,1:n) = 0.0_nag_wp
b(1,2) = 6.0_nag_wp

!     Iterate until coefficients of linearized systems converge.
iter = 0
emax = 1.0_nag_wp
iters: Do While (emax>=coeff_tol)
iter = iter + 1
Write (nout,*)
Write (nout,99999) ' Iteration', iter, ' Chebyshev coefficients are'

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d02tgf(n,m,l,x0,x1,k1,kp,c,ldc,coeff,bdyc,w,lw,iw,liw,ifail)

Write (nout,99998)(c(1:k1,j),j=1,n)
emax = 0.0_nag_wp
Do j = 1, n
Do i = 1, k1
emax = max(emax,abs(c(i,j)-b(i,j)))
End Do
End Do
b(1:k1,1:n) = c(1:k1,1:n)
End Do iters

Deallocate (b)

!     Print solution on uniform mesh.
k = 9
ia1 = 1
Write (nout,*)
Write (nout,99999) 'Solution evaluated at', k, '  equally spaced points'
Write (nout,*)
Write (nout,99997) '      X ', (j,j=1,n)
xinc = (x1-x0)/real(k-1,kind=nag_wp)
x = x0
Do i = 1, k
Do j = 1, n
ifail = 0
Call e02akf(k1,x0,x1,c(1,j),ia1,k1,x,y(j),ifail)
End Do
Write (nout,99996) x, (y(j),j=1,n)
x = x + xinc
End Do

99999 Format (1X,A,I3,A)
99998 Format (1X,9F8.4)
99997 Format (1X,A,2('      Y(',I1,')'))
99996 Format (1X,3F10.4)
End Program d02tgfe
```