C  EXAMPLE FOR THE CALL OF SUBROUTINE KEPLEX
C  TWO BODY PROBLEM WITH ECCENTRICITY 'ECCEN' AND
C  FIXED FREQUENCY 'OM'
C  DUE TO STIEFEL/BETTIS
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Y(4)
      COMMON /PAR/ EPAR
      EXTERNAL FCN
C
C
C  PRINT-PARAMETER KFLAG
C  KFLAG=0 : NO OUTPUT
C  KFLAG=1 : INTEGRATION MONITOR
C  KFLAG>1 : ADDITIONALLY INTERMEDIATE SOLUTION POINTS T,Y(1),...
C
      KFLAG=1
C  NUMBER OF INITIAL VALUES IN Y + NUMBER OF INITIAL VALUES IN DY/DT
      N=4
C  STARTING POINT OF INTEGRATION
      X=0.D0
C  ECCENTRICITY
      ECCEN=0.5D0
C  FIXED FREQUENCY
      EPAR=0.0014D0/3.D0
      OM=0.5D0*DSQRT(1.D0-ECCEN+2.D0*EPAR)
C  INITIAL VALUES Y(T)
      Y(1)=1.D0
      Y(2)=0.D0
C  INITIAL VALUES DY/DT
      Y(3)=0.D0
      Y(4)=0.5D0*DSQRT(1.D0+ECCEN)
C  FINAL POINT OF INTEGRATION
      PER=2.D0*4.4395413186376D0
      XEND= 40.0D0*PER
C  DESIRED ACCURACY
      EPS=1.D-12
C  MAXIMUM PERMITTED STEPSIZE
      HMAX=XEND-X
C  INITIAL STEPSIZE GUESS
      H=1.D-3
C
      WRITE(6,101) X,(Y(I),I=1,N)
C  CALL OF KEPLEX WITH THESE PARAMETERS
      CALL KEPLEX (N,FCN,OM,X,Y,XEND,EPS,HMAX,H,KFLAG)
C
      R=DSQRT(Y(1)*Y(1)+Y(2)*Y(2))
      WRITE(6,201) X,(Y(I),I=1,N),R
101   FORMAT('  KEPLEX: TWO BODY PROBLEM', /,
     $       '  ========================', //,
     $       '  INITIAL VALUES  ',D25.16,/,10X,4D25.16)
201   FORMAT('  SOLUTION AT T=  ',D25.16,/,10X,4D25.16,/,
     $       '  RADIUS          ',D25.16)
      STOP
      END
C
C
      SUBROUTINE FCN (N,T,Z,DZ2)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /PAR/ EPAR
      DIMENSION Z(N),DZ2(N)
C
C  TWO- BODY-PROBLEM
C
      RH=Z(1)*Z(1) + Z(2)*Z(2)
C
      R3=RH*RH*RH
      DZ2(1)=-EPAR*Z(1)/R3
      DZ2(2)=-EPAR*Z(2)/R3
C
      RETURN
      END