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

    MODULE d02mvfe_mod

!      D02MVF 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       :: one = 1.0_nag_wp
       REAL (KIND=nag_wp), PARAMETER       :: two = 2.0_nag_wp
       INTEGER, PARAMETER                  :: iset = 1, neq = 6, nin = 5,      &
                                              nout = 6
       INTEGER, PARAMETER                  :: nrw = 50 + 4*neq
       INTEGER, PARAMETER                  :: nwkjac = neq*(neq+1)
       INTEGER, PARAMETER                  :: sdysav = 8
       INTEGER, PARAMETER                  :: ldysav = neq
    CONTAINS
       SUBROUTINE resid(neq,t,y,ydot,r,ires)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          INTEGER, INTENT (INOUT)             :: ires
          INTEGER, INTENT (IN)                :: neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: r(neq)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq), ydot(neq)
!         .. Executable Statements ..
          IF (ires==-1) THEN
             r(1) = -ydot(1)
             r(2) = -ydot(2)
             r(3) = -ydot(3)
             r(4) = -ydot(4)
             r(5:6) = 0.0E0_nag_wp
          ELSE
             r(1) = y(3) - y(6)*y(1) - ydot(1)
             r(2) = y(4) - y(6)*y(2) - ydot(2)
             r(3) = -y(5)*y(1) - ydot(3)
             r(4) = -y(5)*y(2) - one - ydot(4)
             r(5) = y(1)*y(3) + y(2)*y(4)
             r(6) = y(1)**2 + y(2)**2 - one
          END IF
          RETURN
       END SUBROUTINE resid

       SUBROUTINE jac(neq,t,y,ydot,h,d,p)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: d, h, t
          INTEGER, INTENT (IN)                :: neq
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: p(neq,neq)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(neq), ydot(neq)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: hxd
!         .. Executable Statements ..
          hxd = h*d
          p(1,1) = (one+hxd*y(6))
          p(1,3) = -hxd
          p(1,6) = hxd*y(1)
          p(2,2) = (one+hxd*y(6))
          p(2,4) = -hxd
          p(2,6) = hxd*y(2)
          p(3,1) = hxd*y(5)
          p(3,3) = one
          p(3,5) = hxd*y(1)
          p(4,2) = hxd*y(5)
          p(4,4) = one
          p(4,5) = hxd*y(2)
          p(5,1) = -hxd*y(3)
          p(5,2) = -hxd*y(4)
          p(5,3) = -hxd*y(1)
          p(5,4) = -hxd*y(2)
          p(6,1) = -two*hxd*y(1)
          p(6,2) = -two*hxd*y(2)
          RETURN
       END SUBROUTINE jac
    END MODULE d02mvfe_mod

    PROGRAM d02mvfe

!      D02MVF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02mvf, d02nby, d02ngf, d02nsf, nag_wp, x04abf
       USE d02mvfe_mod, ONLY : iset, jac, ldysav, neq, nin, nout, nrw, nwkjac, &
                               one, resid, sdysav
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: h0, hmax, hmin, t, tcrit, tout
       INTEGER                             :: i, ifail, itask, itol, itrace,   &
                                              maxord, maxstp, mxhnil, outchn
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: atol(:), rtol(:), rwork(:),      &
                                              wkjac(:), y(:), ydot(:), ysav(:,:)
       REAL (KIND=nag_wp)                  :: con(3)
       INTEGER                             :: inform(23)
       LOGICAL                             :: lderiv(2)
!      .. Executable Statements ..
       WRITE (nout,*) 'D02MVF Example Program Results'
!      Skip heading in data file
       READ (nin,*)
!      Allocations based on number of differential equations (neq)

       ALLOCATE (atol(neq),rtol(neq),rwork(nrw),wkjac(nwkjac),y(neq), &
          ydot(neq),ysav(ldysav,sdysav))

!      Read algorithmic parameters
       READ (nin,*) maxord, maxstp, mxhnil
       READ (nin,*) hmin, hmax, h0
       READ (nin,*) rtol(1), atol(1)
       READ (nin,*) itrace, itol
       READ (nin,*) t, tout
       READ (nin,*) y(1:neq)
       READ (nin,*) itask

!      Set initial derivatives and other values
       outchn = nout
       CALL x04abf(iset,outchn)
       con(1:3) = 0.0_nag_wp
       ydot(1) = y(3) - y(6)*y(1)
       ydot(2) = y(4) - y(6)*y(2)
       ydot(3) = -y(5)*y(1)
       ydot(4) = -y(5)*y(2) - one
       ydot(5) = -3.0_nag_wp*y(4)
       ydot(6) = 0.0_nag_wp
       tcrit = tout

!      ifail: behaviour on error exit   
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
       ifail = 0
       CALL d02mvf(neq,sdysav,maxord,con,tcrit,hmin,hmax,h0,maxstp,mxhnil, &
          'AVERAGE-L2',rwork,ifail)

       WRITE (nout,*)
       WRITE (nout,99999) 'Pendulum problem with relative tolerance', rtol(1)
       WRITE (nout,99999) '                  and absolute tolerance', atol(1)
       WRITE (nout,99998) (i,i=1,neq)
       WRITE (nout,99997) t, y(1:neq)

       ifail = 0
       CALL d02nsf(neq,neq,'ANALYTIC',nwkjac,rwork,ifail)

       lderiv(1) = .TRUE.
       lderiv(2) = .TRUE.

       ifail = 0
       CALL d02ngf(neq,ldysav,t,tout,y,ydot,rwork,rtol,atol,itol,inform,resid, &
          ysav,sdysav,jac,wkjac,nwkjac,d02nby,lderiv,itask,itrace,ifail)

       WRITE (nout,99997) t, y(1:neq)

99999  FORMAT (1X,A40,' =',1P,E8.1)
99998  FORMAT (/1X,'  t ',3X,6('      y',I1))
99997  FORMAT (1X,F7.4,2X,6(F8.4))
    END PROGRAM d02mvfe