C  SAMPLE PROGRAM FOR RKF45
C
C
      EXTERNAL ORBIT
      DOUBLE PRECISION T,Y(4),TOUT,RELERR,ABSERR
      DOUBLE PRECISION TFINAL,TPRINT,ECC,ALFA,ALFASQ,WORK(27)
      INTEGER IWORK(5), IFLAG, NEQN
      DOUBLE PRECISION DSQRT
      COMMON ALFASQ
      ECC = 0.25
      ALFA = 3.141592653589/4.0
      ALFASQ = ALFA*ALFA
      NEQN = 4
      T = 0.0
      Y(1) = 1.0 - ECC
      Y(2) = 0.0
      Y(3) = 0.0
      Y(4) = ALFA*DSQRT( (1.0 + ECC)/(1.0 - ECC) )
      RELERR = 1.0D-9
      ABSERR = 0.0
      TFINAL = 12.0
      TPRINT = 0.5
      IFLAG = 1
      TOUT = T
   10 CALL RKF45(ORBIT,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK,IWORK)
      WRITE(1,11) T, Y(1), Y(2)
      GO TO (80,20,30,40,50,60,70,80), IFLAG
   20 TOUT = T + TPRINT
      IF (T .LT. TFINAL) GO TO 10
      STOP
   30 WRITE(1,31) RELERR,ABSERR
      GO TO 10
   40 WRITE(1,41)
      GO TO 10
   50 ABSERR = 1.0D-9
      WRITE(1,31) RELERR,ABSERR
      GO TO 10
   60 RELERR = 10.0*RELERR
      WRITE(1,31) RELERR,ABSERR
      IFLAG = 2
      GO TO 10
   70 WRITE(1,71)
      IFLAG = 2
      GO TO 10
   80 WRITE(1,81)
      STOP
C
   11 FORMAT(F5.1, 2F15.9)
   31 FORMAT(17H TOLERANCES RESET, 2E12.3)
   41 FORMAT(11H MANY STEPS)
   71 FORMAT(12H MUCH OUTPUT)
   81 FORMAT(14H IMPROPER CALL)
      END


      SUBROUTINE ORBIT (T, Y, YP)
      DOUBLE PRECISION T, Y(4), YP(4), R, ALFASQ
      DOUBLE PRECISION DSQRT
      COMMON ALFASQ
      R = Y(1)*Y(1) + Y(2)*Y(2)
      R = R*DSQRT(R)/ALFASQ
      YP(1) = Y(3)
      YP(2) = Y(4)
      YP(3) = -Y(1)/R
      YP(4) = -Y(2)/R
      RETURN
      END
