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

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
!     .. Accessibility Statements ..
Private
Public                           :: fcn
!     .. Parameters ..
Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
Integer, Parameter, Public       :: neq = 2, nin = 5, nout = 6
Integer, Parameter, Public       :: 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
!     neq: number of second-order ordinary differential equations
Allocate (rwork(lrwork),thres(neq),thresp(neq),y(neq),ydp(neq),          &
yinit(neq),yp(neq),ypinit(neq),ypwant(nwant),ywant(nwant))

!     Initial conditions

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 one-step 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) Then
Exit loop2
End If

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
```