NAG Library Manual, Mark 28.4
```!   D02AGF Example Program Text
!   Mark 28.4 Release. NAG Copyright 2022.

Module d02agfe_mod

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

!     iprint:    set iprint = 1 for output at each Newton iteration.
!     nin:       the input channel number
!     nout:      the output channel number

!     For Example 1:
!     n_ex1 : number of differential equations,
!     n1_ex1: number of parameters.

!     For Example 2:
!     n_ex2 : number of differential equations,
!     n1_ex2: number of parameters.

!     .. Use Statements ..
Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Accessibility Statements ..
Private
Public                           :: aux1, aux2, bcaux1, bcaux2, prsol1,  &
prsol2, raaux1, raaux2
!     .. Parameters ..
Integer, Parameter               :: iprint = 0
Integer, Parameter, Public       :: n1_ex1 = 2, n1_ex2 = 3, nin = 5,     &
nout = 6, n_ex1 = 2, n_ex2 = 3
Contains
Subroutine aux1(f,y,x,param)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: x
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: f(*)
Real (Kind=nag_wp), Intent (In) :: param(*), y(*)
!       .. Executable Statements ..
f(1) = y(2)
f(2) = (y(1)**3-y(2))/(2.0_nag_wp*x)
Return
End Subroutine aux1
Subroutine raaux1(x0,x1,r,param)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Out) :: r, x0, x1
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: param(*)
!       .. Executable Statements ..
x0 = 0.1_nag_wp
x1 = 16.0_nag_wp
r = 16.0_nag_wp
Return
End Subroutine raaux1
Subroutine bcaux1(g0,g1,param)

!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: g0(*), g1(*)
Real (Kind=nag_wp), Intent (In) :: param(*)
!       .. Local Scalars ..
Real (Kind=nag_wp)             :: z
!       .. Intrinsic Procedures ..
Intrinsic                      :: sqrt
!       .. Executable Statements ..
z = 0.1_nag_wp
g0(1) = 0.1_nag_wp + param(1)*sqrt(z)*0.1_nag_wp + 0.01_nag_wp*z
g0(2) = param(1)*0.05_nag_wp/sqrt(z) + 0.01_nag_wp
g1(1) = 1.0_nag_wp/6.0_nag_wp
g1(2) = param(2)
Return
End Subroutine bcaux1
Subroutine prsol1(param,res,n1,err)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: res
Integer, Intent (In)           :: n1
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: err(n1), param(n1)
!       .. Local Scalars ..
Integer                        :: i
!       .. Executable Statements ..
If (iprint/=0) Then
Write (nout,99999) 'Current parameters   ', (param(i),i=1,n1)
Write (nout,99998) 'Residuals   ', (err(i),i=1,n1)
Write (nout,99998) 'Sum of residuals squared  ', res
Write (nout,*)
End If
Return

99999   Format (1X,A,6(E14.6,2X))
99998   Format (1X,A,6(E12.4,1X))
End Subroutine prsol1
Subroutine aux2(f,y,x,param)

!       .. Parameters ..
Real (Kind=nag_wp), Parameter  :: eps = 2.0E-5_nag_wp
!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: x
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: f(*)
Real (Kind=nag_wp), Intent (In) :: param(*), y(*)
!       .. Local Scalars ..
Real (Kind=nag_wp)             :: c, s
!       .. Intrinsic Procedures ..
Intrinsic                      :: cos, sin
!       .. Executable Statements ..
c = cos(y(3))
s = sin(y(3))
f(1) = s/c
f(2) = -(param(1)*s+eps*y(2)*y(2))/(y(2)*c)
f(3) = -param(1)/(y(2)*y(2))
Return
End Subroutine aux2
Subroutine raaux2(x0,x1,r,param)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (Out) :: r, x0, x1
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: param(*)
!       .. Executable Statements ..
x0 = 0.0_nag_wp
x1 = param(2)
r = param(2)
Return
End Subroutine raaux2
Subroutine bcaux2(g0,g1,param)

!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (Out) :: g0(*), g1(*)
Real (Kind=nag_wp), Intent (In) :: param(*)
!       .. Executable Statements ..
g0(1) = 0.0E0_nag_wp
g0(2) = 500.0E0_nag_wp
g0(3) = 0.5E0_nag_wp
g1(1) = 0.0E0_nag_wp
g1(2) = 450.0E0_nag_wp
g1(3) = param(3)
Return
End Subroutine bcaux2
Subroutine prsol2(param,res,n1,err)

!       .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: res
Integer, Intent (In)           :: n1
!       .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: err(n1), param(n1)
!       .. Local Scalars ..
Integer                        :: i
!       .. Executable Statements ..
If (iprint/=0) Then
Write (nout,99999) 'Current parameters   ', (param(i),i=1,n1)
Write (nout,99998) 'Residuals   ', (err(i),i=1,n1)
Write (nout,99998) 'Sum of residuals squared  ', res
Write (nout,*)
End If
Return

99999   Format (1X,A,6(E14.6,2X))
99998   Format (1X,A,6(E12.4,1X))
End Subroutine prsol2
End Module d02agfe_mod
Program d02agfe

!     D02AGF Example Main Program

!     .. Use Statements ..
Use d02agfe_mod, Only: nout
Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
Implicit None
!     .. Executable Statements ..
Write (nout,*) 'D02AGF Example Program Results'

Call ex1

Call ex2

Contains
Subroutine ex1

!       .. Use Statements ..
Use d02agfe_mod, Only: aux1, bcaux1, n1_ex1, nin, n_ex1, prsol1,       &
raaux1
Use nag_library, Only: d02agf
!       .. Local Scalars ..
Real (Kind=nag_wp)             :: h, r, soler, x, x1, xx
Integer                        :: i, ifail, j, m1
!       .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: c(:,:), e(:), mat(:,:), param(:),   &
parerr(:), wspac1(:), wspac2(:),     &
wspace(:,:)
Real (Kind=nag_wp)             :: dummy(1,1)
!       .. Intrinsic Procedures ..
Intrinsic                      :: real
!       .. Executable Statements ..
!       Skip heading in data file
!       m1: final solution calculated at m1 points in range including
!       end points.
Allocate (c(m1,n_ex1),e(n_ex1),mat(n_ex1,n_ex1),param(n_ex1),          &
parerr(n_ex1),wspac1(n_ex1),wspac2(n_ex1),wspace(n_ex1,9))
!       h: step size estimate,
!       param: initial parameter estimates,
!       parerr: Newton iteration tolerances,
!       soler: bound on the local error.
e(1:n_ex1) = soler
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Case 1'
Write (nout,*)

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d02agf(h,e,parerr,param,c,n_ex1,n1_ex1,m1,aux1,bcaux1,raaux1,     &
prsol1,mat,dummy,wspace,wspac1,wspac2,ifail)

Write (nout,*) 'Final parameters'
Write (nout,99999)(param(i),i=1,n1_ex1)
Write (nout,*)
Write (nout,*) 'Final solution'
Write (nout,*) 'X-value     Components of solution'
Call raaux1(x,x1,r,param)
h = (x1-x)/real(m1-1,kind=nag_wp)
xx = x
Do i = 1, m1
Write (nout,99998) xx, (c(i,j),j=1,n_ex1)
xx = xx + h
End Do

Return

99999   Format (1X,3E16.6)
99998   Format (1X,F7.2,3E13.4)
End Subroutine ex1
Subroutine ex2

!       .. Use Statements ..
Use d02agfe_mod, Only: aux2, bcaux2, n1_ex2, nin, n_ex2, prsol2,       &
raaux2
Use nag_library, Only: d02agf
!       .. Local Scalars ..
Real (Kind=nag_wp)             :: h, r, soler, x, x1, xx
Integer                        :: i, ifail, j, m1
!       .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: c(:,:), e(:), mat(:,:), param(:),   &
parerr(:), wspac1(:), wspac2(:),     &
wspace(:,:)
Real (Kind=nag_wp)             :: dummy(1,1)
!       .. Executable Statements ..
!       m1: final solution calculated at m1 points in range including
!       end points.
Allocate (c(m1,n_ex2),e(n_ex2),mat(n_ex2,n_ex2),param(n_ex2),          &
parerr(n_ex2),wspac1(n_ex2),wspac2(n_ex2),wspace(n_ex2,9))
Write (nout,*)
Write (nout,*)
Write (nout,*) 'Case 2'
Write (nout,*)
!       h: step size estimate,
!       param: initial parameter estimates,
!       parerr: Newton iteration tolerances,
!       soler: bound on the local error.
e(1:n_ex2) = soler

!       ifail: behaviour on error exit
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call d02agf(h,e,parerr,param,c,n_ex2,n1_ex2,m1,aux2,bcaux2,raaux2,     &
prsol2,mat,dummy,wspace,wspac1,wspac2,ifail)

Write (nout,*) 'Final parameters'
Write (nout,99999)(param(i),i=1,n_ex2)
Write (nout,*)
Write (nout,*) 'Final solution'
Write (nout,*) 'X-value     Components of solution'
Call raaux2(x,x1,r,param)
h = (x1-x)/5.0E0_nag_wp
xx = x
Do i = 1, 6
Write (nout,99998) xx, (c(i,j),j=1,n_ex2)
xx = xx + h
End Do

Return

99999   Format (1X,3E16.6)
99998   Format (1X,F7.0,3E13.4)
End Subroutine ex2
End Program d02agfe
```