```!   D02MVF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

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 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
!     Allocations based on number of differential equations (neq)

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

!     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,  &