!   D02GAF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    MODULE d02gafe_mod

!      Data for D02GAF example program

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       REAL (KIND=nag_wp), PARAMETER       :: one = 1.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: zero = 0.0_nag_wp
       INTEGER, PARAMETER                  :: iset = 1, n = 3, nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: beta
    CONTAINS
       SUBROUTINE fcn(x,y,f)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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) - beta*(1.0E0_nag_wp-y(2)*y(2))
          RETURN
       END SUBROUTINE fcn
    END MODULE d02gafe_mod
    PROGRAM d02gafe

!      D02GAF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02gaf, nag_wp, x04abf
       USE d02gafe_mod, ONLY : beta, fcn, iset, n, nin, nout, one, zero
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: a, b, 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 Functions ..
       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
       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
LOOP:  DO k = 1, 2
          SELECT CASE (k)
          CASE (1)
             beta = zero
          CASE (2)
             beta = 0.2_nag_wp
          END SELECT

!         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) WRITE (nout,99999) 'Problem with BETA = ', beta
          IF (ifail==0 .OR. ifail==3) THEN
             WRITE (nout,*)
             IF (ifail==3) WRITE (nout,*) ' IFAIL = 3'
             WRITE (nout,99998) np
             WRITE (nout,99997)
             WRITE (nout,99996) (x(i),(y(j,i),j=1,n),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