!   D02UEF Example Program Text
!   Mark 29.1 Release. NAG Copyright 2023.

    Module d02uefe_mod

!     D02UEF 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                           :: bndary, exact, pdedef
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: four = 4.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: three = 3.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: two = 2.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: m = 3, nin = 5, nout = 6
      Logical, Parameter, Public       :: reqerr = .False.
!     .. Local Scalars ..
      Real (Kind=nag_wp), Public, Save :: a, b
    Contains
      Function exact(x,q)

!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: exact
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
        Integer, Intent (In)           :: q
!       .. Intrinsic Procedures ..
        Intrinsic                      :: cos, sin
!       .. Executable Statements ..
        Select Case (q)
        Case (0)
          exact = cos(x)
        Case (1)
          exact = -sin(x)
        Case (2)
          exact = -cos(x)
        Case (3)
          exact = sin(x)
        End Select
      End Function exact
      Subroutine bndary(m,y,bmat,bvec)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: m
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: bmat(m,m+1), bvec(m), y(m)
!       .. Executable Statements ..
!       Boundary condition on left side of domain
        y(1:2) = a
        y(3) = b
!       Set up Dirichlet condition using exact solution
        bmat(1:m,1:m+1) = zero
        bmat(1:3,1) = one
        bmat(2:3,2) = two
        bmat(2:3,3) = three
        bvec(1) = zero
        bvec(2) = two
        bvec(3) = -two
        Return
      End Subroutine bndary
      Subroutine pdedef(m,f)

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: m
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(m+1)
!       .. Executable Statements ..
        f(1) = one
        f(2) = two
        f(3) = three
        f(4) = four
        Return
      End Subroutine pdedef
    End Module d02uefe_mod
    Program d02uefe

!     D02UEF Example Main Program

!     .. Use Statements ..
      Use d02uefe_mod, Only: a, b, bndary, exact, m, nin, nout, pdedef,        &
                             reqerr, two, zero
      Use nag_library, Only: d02uaf, d02ubf, d02ucf, d02uef, nag_wp, x01aaf,   &
                             x02ajf
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: pi, resid, teneps
      Integer                          :: i, ifail, n, q, q1
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: bmat(m,m+1), bvec(m), f(m+1),        &
                                          uerr(m+1), y(m)
      Real (Kind=nag_wp), Allocatable  :: c(:), f0(:), u(:,:), uc(:,:), x(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, cos, int, max, sin
!     .. Executable Statements ..
      Write (nout,*) ' D02UEF Example Program Results '
      Write (nout,*)

      Read (nin,*)
      Read (nin,*) n

      Allocate (u(n+1,m+1),f0(n+1),c(n+1),uc(n+1,m+1),x(n+1))

!     Set up domain, boundary conditions and definition
      pi = x01aaf(zero)
      a = -pi/two
      b = pi/two
      Call bndary(m,y,bmat,bvec)
      Call pdedef(m,f)

!     Set up solution grid.
      ifail = 0
      Call d02ucf(n,a,b,x,ifail)

!     Set up problem right hand sides for grid and transform.
      f0(1:n+1) = two*sin(x(1:n+1)) - two*cos(x(1:n+1))
      ifail = 0
      Call d02uaf(n,f0,c,ifail)

!     Solve in coefficient space.
      ifail = 0
      Call d02uef(n,a,b,m,c,bmat,y,bvec,f,uc,resid,ifail)
!     Evaluate solution and derivatives on Chebyshev grid.
      Do q = 0, m
        ifail = 0
        Call d02ubf(n,a,b,q,uc(1,q+1),u(1,q+1),ifail)
      End Do
!     Print solution
      Write (nout,*) ' Numerical Solution U and its first three derivatives'
      Write (nout,*)
      Write (nout,99999)
      Write (nout,99998)(x(i),u(i,1:m+1),i=1,n+1)

      If (reqerr) Then
        uerr(1:m+1) = zero
        Do i = 1, n + 1
          Do q = 0, m
            q1 = q + 1
            uerr(q1) = max(uerr(q1),abs(u(i,q1)-exact(x(i),q)))
          End Do
        End Do
        teneps = 10.0_nag_wp*x02ajf()
        Write (nout,'(//)')
        Write (nout,99997)(q,10*(int(uerr(q+1)/teneps)+1),q=0,m)
      End If

99999 Format (1X,T8,'X',T18,'U',T28,'Ux',T37,'Uxx',T47,'Uxxx')
99998 Format (1X,5F10.4)
99997 Format (1X,'Error in the order ',I1,' derivative of U is < ',I8,         &
        ' * machine precision.')

    End Program d02uefe