    MODULE ivp_ode_rk_ex03_mod

      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Default Accessibility ..
      PUBLIC
      ! .. Intrinsic Functions ..
      INTRINSIC KIND
      ! .. Parameters ..
      INTEGER, PARAMETER :: wp = KIND(1.0D0)

    CONTAINS

      FUNCTION f(t,y)

        ! .. Implicit None Statement ..
        IMPLICIT NONE
        ! .. Intrinsic Functions ..
        INTRINSIC SIZE, SQRT
        ! .. Scalar Arguments ..
        REAL (wp), INTENT (IN) :: t
        ! .. Array Arguments ..
        REAL (wp), INTENT (IN) :: y(:)
        ! .. Function Return Value ..
        REAL (wp) :: f(SIZE(y))
        ! .. Local Scalars ..
        REAL (wp) :: r3
        ! .. Executable Statements ..
        r3 = (SQRT(y(1)**2+y(2)**2))**3
        f(1) = y(3)
        f(2) = y(4)
        f(3) = -y(1)/r3
        f(4) = -y(2)/r3

      END FUNCTION f

    END MODULE ivp_ode_rk_ex03_mod

    PROGRAM nag_ivp_ode_rk_ex03

      ! Example Program Text for nag_ivp_ode_rk
      ! NAG fl90, Release 4. NAG Copyright 2000.

      ! .. Use Statements ..
      USE nag_ivp_ode_rk, ONLY : nag_rk_comm_wp => nag_rk_comm_dp, &
       nag_rk_setup, nag_rk_step, nag_rk_info, nag_rk_reset_end, &
       nag_deallocate
      USE nag_math_constants, ONLY : nag_pi
      USE nag_examples_io, ONLY : nag_std_out
      USE ivp_ode_rk_ex03_mod, ONLY : wp, f
      ! .. Implicit None Statement ..
      IMPLICIT NONE
      ! .. Intrinsic Functions ..
      INTRINSIC SQRT
      ! .. Parameters ..
      INTEGER, PARAMETER :: n = 4
      ! .. Local Scalars ..
      INTEGER :: i, j, method, nout, num_fun_eval
      REAL (wp) :: alpha, pi, tinc, tol, t_end, t_final, t_now, t_start
      LOGICAL :: call_step
      TYPE (nag_rk_comm_wp) :: comm
      ! .. Local Arrays ..
      REAL (wp) :: thresh(n), yp_now(n), y_now(n), y_start(n)
      ! .. Executable Statements ..
      WRITE (nag_std_out,*) 'Example Program Results for nag_ivp_ode_rk_ex03'

      pi = nag_pi(0.0_wp)
      alpha = 0.7_wp

      t_end = 2.0_wp*pi
      t_final = 6.0_wp*pi
      thresh = 1.0E-10_wp
      t_start = 0.0_wp
      y_start = (/ 1.0_wp - alpha, 0.0_wp, 0.0_wp, SQRT((1.0_wp+alpha)/(1.0_wp &
       -alpha)) /)
      method = 3
      call_step = .TRUE.

      nout = 6
      tinc = t_final/nout

      DO i = 1, 2
        IF (i==1) tol = 1.0E-4_wp
        IF (i==2) tol = 4.0E-5_wp
        j = nout - 1
        t_end = t_final - j*tinc

        CALL nag_rk_setup(t_start,y_start,t_end,tol,thresh,comm,method=method, &
         call_step=call_step)

        WRITE (nag_std_out,'(/1X,A,1PE8.1)') 'Calculation with tol = ', tol
        WRITE (nag_std_out,*) &
         '   t         y1         y2         y3         y4'
        WRITE (nag_std_out,'(1X,F6.3,4(4X,F7.4))') t_start, y_start(:)

inner:  DO

          CALL nag_rk_step(f,t_now,y_now,yp_now,comm)

          IF (t_now>=t_end) THEN
            WRITE (nag_std_out,'(1X,F6.3,4(4X,F7.4))') t_now, y_now(:)
            IF (t_now>=t_final) EXIT inner
            j = j - 1
            t_end = t_final - j*tinc

            CALL nag_rk_reset_end(t_end,comm)

          END IF

        END DO inner

        CALL nag_rk_info(comm,num_fun_eval=num_fun_eval)

        WRITE (nag_std_out,'(1X,A,I10)') &
         'The cost of integration in evaluations of f = ', num_fun_eval
      END DO

      CALL nag_deallocate(comm)    ! Free structure allocated by NAG fl90

    END PROGRAM nag_ivp_ode_rk_ex03
