!   D02MVF Example Program Text
!   Mark 28.3 Release. NAG Copyright 2022.

    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
!     .. Accessibility Statements ..
      Private
      Public                           :: jac, resid
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: two = 2.0_nag_wp
      Integer, Parameter, Public       :: iset = 1, neq = 6, nin = 5, nout = 6
      Integer, Parameter, Public       :: nrw = 50 + 4*neq
      Integer, Parameter, Public       :: nwkjac = neq*(neq+1)
      Integer, Parameter, Public       :: sdysav = 8
      Integer, Parameter, Public       :: ldysav = neq
    Contains
      Subroutine resid(neq,t,y,ydot,r,ires)

!       .. 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)

!       .. 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 d02mvfe_mod, Only: iset, jac, ldysav, neq, nin, nout, nrw, nwkjac,   &
                             one, resid, sdysav
      Use nag_library, Only: d02mvf, d02nby, d02ngf, d02nsf, nag_wp, x04abf
!     .. 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