!   D02GAF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.

    Module d02gafe_mod

!     Data for D02GAF example program

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: fcn
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: iset = 1, n = 4, nin = 5, nout = 6
    Contains
      Subroutine fcn(x,y,f)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(*)
        Real (Kind=nag_wp), Intent (In) :: y(*)
!       .. Executable Statements ..
        f(1) = y(2)
        f(2) = y(3)
        f(3) = -y(1)*y(3) - y(4)*(1.0E0_nag_wp-y(2)*y(2))
        f(4) = 0.0_nag_wp
        Return
      End Subroutine fcn
    End Module d02gafe_mod
    Program d02gafe

!     D02GAF Example Main Program

!     .. Use Statements ..
      Use d02gafe_mod, Only: fcn, iset, n, nin, nout, one, zero
      Use nag_library, Only: d02gaf, nag_wp, x04abf
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: a, b, beta, h, tol
      Integer                          :: i, ifail, j, k, liw, lw, mnp, np,    &
                                          outchn
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: u(:,:), v(:,:), w(:), x(:), y(:,:)
      Integer, Allocatable             :: iw(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: real
!     .. Executable Statements ..
      Write (nout,*) 'D02GAF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
!     n: number of differential equations
!     mnp: maximum permitted number of mesh points.
      Read (nin,*) mnp
      liw = mnp*(2*n+1) + n*n + 4*n + 2
      lw = mnp*(3*n*n+6*n+2) + 4*n*n + 4*n
      Allocate (iw(liw),u(n,2),v(n,2),w(lw),x(mnp),y(n,mnp))
!     tol: positive absolute error tolerance
!     np : determines whether a default or user-supplied mesh is used.
!     a  : left-hand boundary point, b: right-hand boundary point.
      Read (nin,*) tol
      Read (nin,*) np
      Read (nin,*) a, b
      outchn = nout
      Call x04abf(iset,outchn)
      beta = zero
      u(1:n,1:2) = zero
      v(1:n,1:2) = zero
      v(1,2) = one
      v(3,1) = one
      v(3,2) = one
      v(4,2) = one
      u(2,2) = one
      u(1,2) = b
      x(1) = a
      h = (b-a)/real(np-1,kind=nag_wp)
      Do i = 2, np - 1
        x(i) = x(i-1) + h
      End Do
      x(np) = b
      beta = zero
loop: Do k = 1, 2
        u(4,1) = beta
        u(4,2) = beta
!       ifail: behaviour on error exit
!              =1 for quiet-soft exit
!       * Set ifail to 111 to obtain monitoring information *
        ifail = 1
        Call d02gaf(u,v,n,a,b,tol,fcn,mnp,x,y,np,w,lw,iw,liw,ifail)

        If (ifail>=0) Then
          Write (nout,99999) 'Problem with BETA = ', beta
        End If
        If (ifail==0 .Or. ifail==3) Then
          Write (nout,*)
          If (ifail==3) Then
            Write (nout,*) ' IFAIL = 3'
          End If
          Write (nout,99998) np
          Write (nout,99997)
          Write (nout,99996)(x(i),(y(j,i),j=1,3),i=1,np)
          beta = beta + 0.2E0_nag_wp
        Else
          Write (nout,99995) ifail
          Exit loop
        End If
      End Do loop

99999 Format (/,1X,A,F7.2)
99998 Format (1X,'Solution on final mesh of ',I2,' points')
99997 Format (1X,'      X(I)        Y1(I)        Y2(I)        Y3(I)')
99996 Format (1X,F11.3,3F13.4)
99995 Format (1X,/,1X,' ** D02GAF returned with IFAIL = ',I5)
    End Program d02gafe