!   D02LAF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    Module d02lafe_mod

!     D02LAF 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        :: zero = 0.0_nag_wp
      Integer, Parameter                   :: neq = 2, nin = 5, nout = 6
      Integer, Parameter                   :: lrwork = 16 + 20*neq
    Contains
      Subroutine fcn(neq,t,y,f)

!       Derivatives for two body problem in  y'' = f(t,y)  form

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: t
        Integer, Intent (In)                 :: neq
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: f(neq)
        Real (Kind=nag_wp), Intent (In)      :: y(neq)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: r
!       .. Intrinsic Procedures ..
        Intrinsic                            :: sqrt
!       .. Executable Statements ..
        r = sqrt(y(1)**2+y(2)**2)**3
        f(1) = -y(1)/r
        f(2) = -y(2)/r
        Return
      End Subroutine fcn
    End Module d02lafe_mod

    Program d02lafe

!     D02LAF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: d02laf, d02lxf, d02lyf, d02lzf, nag_wp
      Use d02lafe_mod, Only: fcn, lrwork, neq, nin, nout, zero
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: h, hnext, hstart, hused, t,      &
                                              tend, tinc, tnext, tol, tstart
      Integer                              :: i, ifail, itol, maxstp, natt,    &
                                              nfail, nsucc, nwant
      Logical                              :: high, onestp, start
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: rwork(:), thres(:), thresp(:),   &
                                              y(:), ydp(:), yinit(:), yp(:),   &
                                              ypinit(:), ypwant(:), ywant(:)
!     .. Executable Statements ..
      Write (nout,*) 'D02LAF Example Program Results'
!     Skip heading in data file
      Read (nin,*)
!     neq: number of second-order ordinary differential equations
      Read (nin,*) nwant
      Allocate (rwork(lrwork),thres(neq),thresp(neq),y(neq),ydp(neq), &
        yinit(neq),yp(neq),ypinit(neq),ypwant(nwant),ywant(nwant))
      Read (nin,*) high, onestp
      Read (nin,*) tinc

!     Initial conditions
      Read (nin,*) tstart, tend
      Read (nin,*) yinit(1:neq)
      Read (nin,*) ypinit(1:neq)

loop1: Do itol = 4, 5
        tol = 10.0_nag_wp**(-itol)
        Write (nout,*)

!       Call D02LXF with default THRES,THRESP,MAXSTP and H

        thres(1) = zero
        thresp(1) = zero
        h = zero
        maxstp = 0
        start = .True.

!       ifail: behaviour on error exit   
!              =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
        ifail = 0
        Call d02lxf(neq,h,tol,thres,thresp,maxstp,start,onestp,high,rwork, &
          lrwork,ifail)

        Write (nout,99999) 'Calculation with TOL = ', tol
        Write (nout,99995)(i,i=1,neq)

!       Set initial values

        y(1:neq) = yinit(1:neq)
        yp(1:neq) = ypinit(1:neq)
        t = tstart
        tnext = t + tinc
        Write (nout,99998) t, y(1:neq)

!       Loop point for onestep mode
loop2:  Do

          ifail = -1
          Call d02laf(fcn,neq,t,tend,y,yp,ydp,rwork,lrwork,ifail)

          If (ifail>0) Then
            Write (nout,99997) ifail, t
            Exit loop1
          End If

!         Loop point for interpolation
          Do While (tnext<=t)

            ifail = 0
            Call d02lzf(neq,t,y,yp,neq,tnext,ywant,ypwant,rwork,lrwork,ifail)

            Write (nout,99998) tnext, ywant(1:neq)
            tnext = tnext + tinc
          End Do

          If (t>=tend) Exit loop2

        End Do loop2

        ifail = 0
        Call d02lyf(neq,hnext,hused,hstart,nsucc,nfail,natt,thres,thresp, &
          rwork,lrwork,ifail)

        Write (nout,*)
        Write (nout,99996) ' Number of successful steps = ', nsucc
        Write (nout,99996) ' Number of   failed   steps = ', nfail
      End Do loop1

99999 Format (1X,A,1P,E9.1)
99998 Format (1X,F5.1,2(2X,F9.5))
99997 Format (/1X,'D02LAF returned with IFAIL = ',I2,'  at T = ',1P,E10.3)
99996 Format (1X,A,I5)
99995 Format (/'   T ',2('     Y(',I1,')  '))
    End Program d02lafe
