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

    Module d02txfe_mod

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

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter        :: half = 0.5_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        :: two = 2.0_nag_wp
      Real (Kind=nag_wp), Parameter        :: xsplit = 30.0_nag_wp
      Integer, Parameter                   :: m1 = 3, m2 = 2, mmax = 3,        &
                                              neq = 2, nin = 5, nlbc = 3,      &
                                              nleft = 15, nout = 6, nrbc = 2,  &
                                              nright = 10
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: el, en, s
    Contains
      Subroutine ffun(x,y,neq,m,f)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: x
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: f(neq)
        Real (Kind=nag_wp), Intent (In)      :: y(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: t1, y11, y20
!       .. Executable Statements ..
        t1 = half*(three-en)*y(1,0)
        y11 = y(1,1)
        y20 = y(2,0)
        f(1) = (el**3)*(one-y20**2) + (el**2)*s*y11 - el*(t1*y(1,2)+en*y11**2)
        f(2) = (el**2)*s*(y20-one) - el*(t1*y(2,1)+(en-one)*y11*y20)
        Return
      End Subroutine ffun
      Subroutine fjac(x,y,neq,m,dfdy)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: x
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: dfdy(neq,neq,0:*)
        Real (Kind=nag_wp), Intent (In)      :: y(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Executable Statements ..
        dfdy(1,2,0) = -two*el**3*y(2,0)
        dfdy(1,1,0) = -el*half*(three-en)*y(1,2)
        dfdy(1,1,1) = el**2*s - el*two*en*y(1,1)
        dfdy(1,1,2) = -el*half*(three-en)*y(1,0)
        dfdy(2,2,0) = el**2*s - el*(en-one)*y(1,1)
        dfdy(2,2,1) = -el*half*(three-en)*y(1,0)
        dfdy(2,1,0) = -el*half*(three-en)*y(2,1)
        dfdy(2,1,1) = -el*(en-one)*y(2,0)
        Return
      End Subroutine fjac
      Subroutine gafun(ya,neq,m,nlbc,ga)

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: neq, nlbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: ga(nlbc)
        Real (Kind=nag_wp), Intent (In)      :: ya(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Executable Statements ..
        ga(1) = ya(1,0)
        ga(2) = ya(1,1)
        ga(3) = ya(2,0)
        Return
      End Subroutine gafun
      Subroutine gbfun(yb,neq,m,nrbc,gb)

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: neq, nrbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: gb(nrbc)
        Real (Kind=nag_wp), Intent (In)      :: yb(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Executable Statements ..
        gb(1) = yb(1,1)
        gb(2) = yb(2,0) - one
        Return
      End Subroutine gbfun
      Subroutine gajac(ya,neq,m,nlbc,dgady)

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: neq, nlbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: dgady(nlbc,neq,0:*)
        Real (Kind=nag_wp), Intent (In)      :: ya(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Executable Statements ..
        dgady(1,1,0) = one
        dgady(2,1,1) = one
        dgady(3,2,0) = one
        Return
      End Subroutine gajac
      Subroutine gbjac(yb,neq,m,nrbc,dgbdy)

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: neq, nrbc
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: dgbdy(nrbc,neq,0:*)
        Real (Kind=nag_wp), Intent (In)      :: yb(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Executable Statements ..
        dgbdy(1,1,1) = one
        dgbdy(2,2,0) = one
        Return
      End Subroutine gbjac
      Subroutine guess(x,neq,m,y,dym)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: x
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: dym(neq)
        Real (Kind=nag_wp), Intent (Inout)   :: y(neq,0:*)
        Integer, Intent (In)                 :: m(neq)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: ex, expmx
!       .. Intrinsic Procedures ..
        Intrinsic                            :: exp
!       .. Executable Statements ..
        ex = x*el
        expmx = exp(-ex)
        y(1,0) = -ex**2*expmx
        y(1,1) = (-two*ex+ex**2)*expmx
        y(1,2) = (-two+4.0E0_nag_wp*ex-ex**2)*expmx
        y(2,0) = one - expmx
        y(2,1) = expmx
        dym(1) = (6.0E0_nag_wp-6.0E0_nag_wp*ex+ex**2)*expmx
        dym(2) = -expmx
        Return
      End Subroutine guess
    End Module d02txfe_mod
    Program d02txfe

!     D02TXF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02tkf, d02tvf, d02txf, d02tyf, d02tzf, nag_wp
      Use d02txfe_mod, Only: el, en, ffun, fjac, gafun, gajac, gbfun, gbjac,   &
                             guess, m1, m2, mmax, neq, nin, nlbc, nleft, nout, &
                             nrbc, nright, one, s, two, xsplit
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: dx, el_init, ermx, s_init, xx
      Integer                              :: i, iermx, ifail, ijermx, j,      &
                                              liwork, lrwork, mxmesh, ncol,    &
                                              ncont, nmesh
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: mesh(:), rwork(:)
      Real (Kind=nag_wp)                   :: tol(neq), y(neq,0:mmax-1)
      Integer, Allocatable                 :: ipmesh(:), iwork(:)
      Integer                              :: m(neq)
!     .. Intrinsic Procedures ..
      Intrinsic                            :: min, real
!     .. Executable Statements ..
      Write (nout,*) 'D02TXF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
!     Read method parameters
      Read (nin,*) ncol, nmesh, mxmesh
      Read (nin,*) tol(1:neq)

      liwork = mxmesh*(11*neq+6)
      lrwork = mxmesh*(109*neq**2+78*neq+7)
      Allocate (mesh(mxmesh),rwork(lrwork),ipmesh(mxmesh),iwork(liwork))
!     Read problem (initial) parameters
      Read (nin,*) en, el_init, s_init
!     Initialise data
      el = el_init
      s = s_init
      m(1) = m1
      m(2) = m2

      dx = one/real(nmesh-1,kind=nag_wp)
      mesh(1) = 0.0_nag_wp
      Do i = 2, nmesh - 1
        mesh(i) = mesh(i-1) + dx
      End Do
      mesh(nmesh) = 1.0_nag_wp

      ipmesh(1) = 1
      ipmesh(2:nmesh-1) = 2
      ipmesh(nmesh) = 1

!     ifail: behaviour on error exit   
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
      ifail = 0
      Call d02tvf(neq,m,nlbc,nrbc,ncol,tol,mxmesh,nmesh,mesh,ipmesh,rwork, &
        lrwork,iwork,liwork,ifail)

!     Initialize number of continuation steps in el and s
      Read (nin,*) ncont
cont: Do j = 1, ncont
        Write (nout,99997) tol(1), el, s

!       Solve
        ifail = -1
        Call d02tkf(ffun,fjac,gafun,gbfun,gajac,gbjac,guess,rwork,iwork,ifail)
        If (ifail/=0) Exit cont

!       Extract mesh
        ifail = 0
        Call d02tzf(mxmesh,nmesh,mesh,ipmesh,ermx,iermx,ijermx,rwork,iwork, &
          ifail)

        Write (nout,99996) nmesh, ermx, iermx, ijermx
!       Print solution components on mesh
        Write (nout,99999)

!       Left side domain [0,xsplit], evaluate at nleft+1 uniform grid points.
        dx = xsplit/real(nleft,kind=nag_wp)/el
        xx = 0.0_nag_wp
        Do i = 0, nleft
          ifail = 0
          Call d02tyf(xx,y,neq,mmax,rwork,iwork,ifail)
          Write (nout,99998) xx*el, y(1,0), y(2,0)
          xx = xx + dx
        End Do

!       Right side domain (xsplit,L], evaluate at nright uniform grid points.
        dx = (el-xsplit)/real(nright,kind=nag_wp)/el
        xx = xsplit/el
        Do i = 1, nright
          xx = min(1.0_nag_wp,xx+dx)
          ifail = 0
          Call d02tyf(xx,y,neq,mmax,rwork,iwork,ifail)
          Write (nout,99998) xx*el, y(1,0), y(2,0)
        End Do

!       Select mesh for continuation and update continuation parameters.
        If (j<ncont) Then
          el = two*el
          s = 0.6_nag_wp*s
          nmesh = (nmesh+1)/2
          ifail = 0
          Call d02txf(mxmesh,nmesh,mesh,ipmesh,rwork,iwork,ifail)
        End If
      End Do cont

99999 Format (/' Solution on original interval:'/6X,'x',8X,'f',10X,'g')
99998 Format (1X,F8.2,2(1X,F10.4))
99997 Format (//' Tolerance = ',E8.1,'  L = ',F8.3,'  S =',F7.4)
99996 Format (/' Used a mesh of ',I4,' points'/' Maximum error = ',E10.2, &
        '  in interval ',I4,' for component ',I4)
    End Program d02txfe
