!   D02JBF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    Module d02jbfe_mod

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

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter                   :: nin = 5, nout = 6
    Contains
      Function cf(i,j,x)

!       .. Function Return Value ..
        Real (Kind=nag_wp)                   :: cf
!       .. Parameters ..
        Integer, Parameter                   :: n = 2
        Real (Kind=nag_wp), Parameter        ::  a(n,n) =                      &
      reshape((/0.0E0_nag_wp,-1.0E0_nag_wp,1.0E0_nag_wp,0.0E0_nag_wp/),(/n,n/))
        Real (Kind=nag_wp), Parameter        ::                                &
                                           r(n) = (/0.0E0_nag_wp,1.0E0_nag_wp/)
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: x
        Integer, Intent (In)                 :: i, j
!       .. Intrinsic Procedures ..
        Intrinsic                            :: reshape
!       .. Executable Statements ..
        If (j>0) cf = a(i,j)
        If (j==0) cf = r(i)
        Return
      End Function cf

      Subroutine bc(i,j,rhs)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: rhs
        Integer, Intent (In)                 :: i
        Integer, Intent (Out)                :: j
!       .. Executable Statements ..
        rhs = 0.0E0_nag_wp
        If (i>1) Then
          j = -1
        Else
          j = 1
        End If
        Return
      End Subroutine bc
    End Module d02jbfe_mod

    Program d02jbfe

!     D02JBF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02jbf, e02akf, nag_wp
      Use d02jbfe_mod, Only: bc, cf, nin, nout
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: dx, x, x0, x1
      Integer                              :: i, ia1, ifail, j, k1, k1max, kp, &
                                              kpmax, ldc, liw, lw, m, n
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: c(:,:), w(:), y(:)
      Integer, Allocatable                 :: iw(:)
!     .. Intrinsic Procedures ..
      Intrinsic                            :: real
!     .. Executable Statements ..
      Write (nout,*) 'D02JBF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
!     n: order of the system of differential equations
!     k1: number of coefficients to be returned
!     kp: number of collocation points
      Read (nin,*) n, k1max, kpmax
      ldc = k1max
      liw = n*(k1max+2)
      lw = 2*n*(kpmax+1)*(n*k1max+1) + 7*n*k1max
      Allocate (iw(liw),c(ldc,n),w(lw),y(n))
!     x0: left-hand boundary, x1: right-hand boundary.
      Read (nin,*) x0, x1
      Write (nout,*)
      Write (nout,*) ' KP  K1   Chebyshev coefficients'
      Do kp = 10, kpmax, 5
        Do k1 = 4, k1max, 2

!         ifail: behaviour on error exit   
!                =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
          ifail = 0
          Call d02jbf(n,cf,bc,x0,x1,k1,kp,c,ldc,w,lw,iw,liw,ifail)

          Write (nout,99999) kp, k1, c(1:k1,1)
          Write (nout,99998)(c(1:k1,j),j=2,n)
          Write (nout,*)
        End Do
      End Do
      k1 = 8
      m = 9
      ia1 = 1
      Write (nout,99997) 'Last computed solution evaluated at', m, &
        '  equally spaced points'
      Write (nout,*)
      Write (nout,99996) '      X ', (j,j=1,n)
      dx = (x1-x0)/real(m-1,kind=nag_wp)
      x = x0
      Do i = 1, m
        Do j = 1, n

          ifail = 0
          Call e02akf(k1,x0,x1,c(1,j),ia1,ldc,x,y(j),ifail)

        End Do
        Write (nout,99995) x, y(1:n)
        x = x + dx
      End Do

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