      SUBROUTINE KEPLEX (N,FCN,OM,T,Y,TEND,EPS,HMAX,H,KFLAG)
C
C* Begin Prologue KEPLEX
C
C  ---------------------------------------------------------------------
C
C* Title
C
C    Explicit extrapolation integrator for non-stiff second-order
C    systems of ordinary differential equations of the form
C
C             y" + (om**2) * y  = f(t,y) ,
C                                          om : given fixed frequency
C
C* Written by        P. Deuflhard, U. Nowak, U. Poehle
C* Purpose           Solution of systems of initial value problems
C* Method            Kepler discretization with h**2-extrapolation
C* Category          i1a1c2. - System of nonstiff second order
C                              differential equations
C* Keywords          extrapolation, ODE, Kepler discretization,
C                    nonstiff
C* Version           1.0 , August 1982
C* Latest Change     February 1991
C* Library           CodeLib
C* Code              Fortran 77
C                    Double Precision
C* Environment       Standard version for FORTRAN77 environments on
C                    PCs, workstations, and hosts
C* Copyright     (c) Konrad-Zuse-Zentrum fuer Informationstechnik
C                    Berlin (ZIB)
C                    Takustrasse 7, D-14195 Berlin-Dahlem
C                    phone : + 49/30/84185-0
C                    fax   : + 49/30/84185-125
C* Contact           Rainald Ehrig
C                    ZIB, Numerical Analysis and Modelling
C                    phone : + 49/30/84185-282
C                    fax   : + 49/30/84185-107
C                    e-mail: ehrig@zib.de
C
C  ---------------------------------------------------------------------
C
C* Licence
C  -------
C
C  You may use or modify this code for your own non-commercial
C  purposes for an unlimited time. 
C  In any case you should not deliver this code without a special 
C  permission of ZIB.
C  In case you intend to use the code commercially, we oblige you
C  to sign an according licence agreement with ZIB.
C
C
C* Warranty
C  --------
C 
C  This code has been tested up to a certain level. Defects and
C  weaknesses, which may be included in the code, do not establish
C  any warranties by ZIB. ZIB does not take over any liabilities
C  which may follow from aquisition or application of this code.
C
C
C* Software status 
C  ---------------
C
C  This code is under care of ZIB and belongs to ZIB software
C  class III.
C
C
C  ---------------------------------------------------------------------
C
C
C  REFERENCES:
C /1/ P. DEUFLHARD:
C     A STUDY OF EXTRAPOLATION METHODS BASED ON MULTISTEP SCHEMES
C     WITHOUT PARASITIC SOLUTIONS
C     J.APPL.MATH.PHYS. (ZAMP) 30, 177-189 (1979)
C /2/ P. DEUFLHARD:
C     KEPLER DISCRETIZATION IN REGULAR CELESTIAL MECHANICS
C     CEL.MECH. 21, 213-223 (1980)
C /3/ P. DEUFLHARD:
C     ORDER AND STEPSIZE CONTROL IN EXTRAPOLATION METHODS
C     TO APPEAR IN NUM.MATH.  (1983)
C
C
C  EXTERNAL SUBROUTINE (TO BE SUPPLIED BY THE USER)
C
C    FCN(NF,T,Y,DY2)    RIGHT-HAND SIDE OF SECOND-ORDER
C                       DIFFERENTIAL EQUATIONS
C      NF (.LE.6)       NUMBER OF SECOND-ORDER ODE'S
C      T                ACTUAL POSITION
C      Y(NF)            VALUES AT T
C      DY2(NF)          SECOND DERIVATIVES AT T
C
C
C  INPUT PARAMETERS (* MARKS TRANSIENT PARAMETERS)
C
C    N=NF+NF (.LE.12)   NUMBER OF SECOND-ORDER DIFFERENTIAL EQUATIONS
C    OM                 FIXED FREQUENCY
C  * T                  STARTING POINT OF INTEGRATION
C  * Y(N)               INITIAL VALUES Y(1),...,Y(NF)
C                       INITIAL DERIVATIVES Y(NF+1),...,Y(N)
C    TEND               PRESCRIBED FINAL POINT OF INTEGRATION
C    EPS                PRESCRIBED RELATIVE PRECISION (.GT.0)
C    HMAX               MAXIMUM PERMITTED STEPSIZE
C  * H                  INITIAL STEPSIZE GUESS
C  * KFLAG              PRINT PARAMETER
C                        0   NO OUTPUT
C                        1   INTEGRATION MONITOR
C                        2   ADDITIONALLY INTERMEDIATE SOLUTION POINTS
C                            T,Y(I),I=1,N
C
C                    (OUTPUT IS WRITTEN ON LOGICAL UNIT LOUT=6)
C
C  OUTPUT PARAMETERS
C
C    T                  ACHIEVED FINAL POINT OF INTEGRATION
C    Y(N)               FINAL VALUES AND DERIVATIVES
C    H                  STEPSIZE PROPOSAL FOR NEXT INTEGRATION STEP
C                       (H.EQ.0. ,IF KEPLEX FAILS TO PROCEED)
C    KFLAG       .GE. 0:SUCCESSFUL INTEGRATION
C                       (KFLAG NOT ALTERED INTERNALLY)
C                .EQ.-1:MORE THAN NSTMAX BASIC INTEGRATION STEPS PER
C                       INTERVAL HAVE BEEN PERFORMED
C                .EQ.-2:MORE THAN JRMAX STEPSIZE REDUCTIONS
C                       OCCURRED PER BASIC INTEGRATION STEP
C                .EQ.-3:STEPSIZE PROPOSAL FOR NEXT BASIC INTEGRATION
C                       TOO SMALL
C                .EQ.-4:N ODD
C
C
C* End Prologue
C  ------------
C
C
C  COMMON /COUNT/ NSTEP,NFCN            (INTERNALLY INITIALIZED)
C
C    NSTEP              NUMBER OF INTEGRATION STEPS
C    NFCN               NUMBER OF FCN-EVALUATIONS
C
C
C
C  CONVERSION TO SINGLE PRECISION:
C  JUST CANCEL THE CSP IN THE CSP COMMENT CARDS AND CONVERT
C  THE ASSOCIATED PRECEDING CARDS TO CDP COMMENT CARDS
C
C
      INTEGER NJ(13),INCR(13),NRED(12)
C     INTEGER NJ(JM),INCR(JM),NRED(KM)
      DOUBLE PRECISION  Y(N),US(12),DUS( 6),DLU( 6),A(13),AL(13,13)
CSP       REAL          Y(N),US(12),DUS( 6),DLU( 6),A(13),AL(13,13)
C                       Y(N),US( N),DUS(NF),DLU(NF),A(JM),AL(JM,JM)
      DOUBLE PRECISION  DLV( 6),UQ( 6),DUQ( 6),DDU( 6)
CSP       REAL          DLV( 6),UQ( 6),DUQ( 6),DDU( 6)
C                       DLV(NF),UQ(NF),DUQ(NF),DDU(NF)
      DOUBLE PRECISION UH( 6),ZH( 6),S(12),DT(12,13),D(13,13)
CSP       REAL         UH( 6),ZH( 6),S(12),DT(12,13),D(13,13)
C                      UH(NF),ZH(NF),S( N),DT( N,JM),D(JM,JM)
      DOUBLE PRECISION B,B1,C,EPH,EPMACH,EPS,ERR,FC,FCM,FCO,FJ,FJ1
CSP       REAL         B,B1,C,EPH,EPMACH,EPS,ERR,FC,FCM,FCO,FJ,FJ1
      DOUBLE PRECISION FMIN,FN,G,GH,GH0,GH1,H,HALF,HMAX,HMAXU
CSP       REAL         FMIN,FN,G,GH,GH0,GH1,H,HALF,HMAX,HMAXU
      DOUBLE PRECISION HR,H1,OM,OMJ,OMJO,ONE,ONE1,Q,RED,RO,SAFE
CSP       REAL         HR,H1,OM,OMJ,OMJO,ONE,ONE1,Q,RED,RO,SAFE
      DOUBLE PRECISION SH,SH1,SMALL,SUM1,SUM2,TA,TEN,TWO,U
CSP       REAL         SH,SH1,SMALL,SUM1,SUM2,TA,TEN,TWO,U
      DOUBLE PRECISION V,W,T,TEND,TEPS,TN,Z,ZERO
CSP       REAL         V,W,T,TEND,TEPS,TN,Z,ZERO
C
      EXTERNAL FCN
C
      COMMON /COUNT/ NSTEP,NFCN
C
      DATA LOUT/6/
C
      DATA ZERO/0.D0/,FMIN/1.D-2/,RO/0.25D0/,HALF/0.5D0/
CSP       DATA ZERO/0.E0/,FMIN/1.E-2/,RO/0.25E0/,HALF/0.5E0/
      DATA ONE/1.D0/,ONE1/1.01D0/,TWO/2.D0/,TEN/1.D1/,SAFE/0.5D0/
CSP       DATA ONE/1.E0/,ONE1/1.01E0/,TWO/2.E0/,TEN/1.E1/,SAFE/0.5E0/
      DATA DT/156*0.D0/
CSP       DATA DT/156*0.E0/
C     DATA DT/N*JM*0.D0/
C
C  STEPSIZE SEQUENCE HA (HARMONIC SEQUENCE)
      DATA NJ/1,2,3,4,5,6,7,8,9,10,11,12,13/
CC
C  MAXIMUM COLUMN NUMBER
      KM=12
C
C  ASSOCIATED MAXIMUM ROW NUMBER
      JM=KM+1
C
C  INTERNAL PARAMETERS
C  STANDARD VALUES FIXED BELOW
C
C
C  MAXIMUM PERMITTED NUMBER OF INTEGRATION STEPS PER INTERVAL
      NSTMAX=10000
C
C  MAXIMUM PERMITTED NUMBER OF STEPSIZE REDUCTIONS
      JRMAX=5
C
C  INITIAL PREPARATIONS
      CALL ZIBCONST(EPMACH,SMALL)
      NSTEP=0
      NFCN=0
      HMAX=DABS(HMAX)
CSP       HMAX=ABS(HMAX)
      HMAXU=HMAX
      NF=N/2
      IF(N.NE.NF+NF) GOTO 54
      FN=DFLOAT(N)
CSP       FN=FLOAT(N)
      TEPS=(DABS(T)+DABS(TEND))*EPMACH
CSP       TEPS=(ABS(T)+ABS(TEND))*EPMACH
      H1=TEND-T
      HR=H
      IF(DABS(H1).LE.TEPS) GOTO 50
CSP       IF(ABS(H1).LE.TEPS) GOTO 50
      Q=H1/H
      IF(Q.LE.EPMACH) GOTO 50
      EPH=RO*EPS
      OMJO=ZERO
      U1=ZERO
      FJ1=DFLOAT(NJ(1))
CSP       FJ1=FLOAT(NJ(1))
      A(1)=FJ1+ONE
      DO 12 J=2,JM
      J1=J-1
      INCR(J1)=0
      NRED(J1)=0
      FJ=DFLOAT(NJ(J))
CSP       FJ=FLOAT(NJ(J))
      V=A(J1)+FJ
      A(J)=V
      DO 10 K=1,J1
      W=FJ/DFLOAT(NJ(K))
CSP       W=FJ/FLOAT(NJ(K))
10    D(J,K)=W*W
      IF(J.EQ.2) GOTO 12
      W=V-FJ1
      DO 11 K1=2,J1
      K=K1-1
      U=(A(K1)-V)/(W*DFLOAT(K+K1))
CSP       U=(A(K1)-V)/(W*FLOAT(K+K1))
      U=EPH**U
      AL(J1,K)=U
      U=U*A(K1)
      IF ( K1.GT.2 .AND. U*ONE1.GT.U1 ) GOTO 13
11    U1=U
12    CONTINUE
      KOH=KM
      JOH=JM
      GOTO 14
13    JOH=MAX0(2,J1-1)
      KOH=JOH-1
14    KM=KOH
      JM=KM+1
      INCR(JM)=-1
      IF(KFLAG.GE.0) WRITE(LOUT,1001) EPS,KM
      IF(KFLAG.GE.1) WRITE(LOUT,1005)
      EPMACH=EPMACH*TEN
15    IF(Q.GE.ONE1) GOTO 16
      H=H1
16    JRED=0
      IF(KFLAG.GT.0) WRITE(LOUT,1009) NSTEP,NFCN,T,K,KOH
      IF(KFLAG.GT.1) WRITE(LOUT,1000) NSTEP,NFCN,T,(Y(I),I=1,N)
      DO 17 K=1,KM
17    INCR(K)=INCR(K)+1
      HMAX=DABS(H1)
CSP       HMAX=ABS(H1)
      IF(HMAXU.LT.HMAX) HMAX=HMAXU
      CALL FCN(NF,T,Y,DLV)
      NFCN=NFCN+1
20    TN=T+H
      IF(H1.EQ.H) TN=TEND
      FCM=DABS(H)/HMAX
CSP       FCM=ABS(H)/HMAX
      IF(FCM.LT.FMIN) FCM=FMIN
C
C  KEPLER DISCRETIZATION
      DO 35 J=1,JM
      M=NJ(J)
      B=H/DFLOAT(M)
CSP       B=H/FLOAT(M)
      G=B*HALF
      IF(OM.NE.ZERO) GOTO 21
      SH1=ZERO
      GH1=B
      GOTO 22
21    CONTINUE
      GH1=DSIN(OM*B)/OM
CSP       GH1=SIN(OM*B)/OM
      SH1=DSIN(OM*G)
CSP       SH1=SIN(OM*G)
      SH1=-TWO*SH1*SH1
22    GH0=GH1*G
      DO 23 I=1,NF
      DUQ(I)=SH1*Y(I)+GH1*Y(NF+I)
      UQ(I)=Y(I)+DUQ(I)
      DDU(I)=GH0*DLV(I)
      DUS(I)=DDU(I)
23    US(I)=UQ(I)+DUS(I)
      SH=SH1*TWO
      GH=GH0*TWO
      M=M-1
      IF(M.EQ.0) GOTO 25
      DO 24 K=1,M
      CALL  FCN(NF,T+B*DFLOAT(K),US,DLU)
CSP       CALL  FCN(NF,T+B*FLOAT(K),US,DLU)
      NFCN=NFCN+1
      DO 24 I=1,NF
      DUQ(I)=DUQ(I)+SH*UQ(I)
      UQ(I)=UQ(I)+DUQ(I)
      DDU(I)=DDU(I)+SH*DUS(I)+GH*DLU(I)
      DUS(I)=DUS(I)+DDU(I)
24    US(I)=UQ(I)+DUS(I)
25    CALL FCN(NF,TN,US,DLU)
      NFCN=NFCN+1
      SUM1=ZERO
      SUM2=ZERO
      DO 27 I=1,NF
      Z=US(I)
      SUM1=SUM1+Z*Z
      US(I)=DUS(I)
      IF(J.GT.1) GOTO 26
      UH(I)=UQ(I)
      ZH(I)=(DUQ(I)+SH1*UQ(I))/GH1
26    US(NF+I)=(DDU(I)+GH0*DLU(I)+SH1*DUS(I))/GH1
      Z=US(NF+I)+ZH(I)
27    SUM2=SUM2+Z*Z
C  SCALING
      SUM1=DSQRT(SUM1)
CSP       SUM1=SQRT(SUM1)
      SUM2=DSQRT(SUM2)
CSP       SUM2=SQRT(SUM2)
      DO 28 I=1,NF
      S(I)=SUM1
28    S(NF+I)=SUM2
C
C  EXTRAPOLATION
C
      ERR=ZERO
      DO 31 I=1,N
      V=DT(I,1)
      C=US(I)
      DT(I,1)=C
      IF(J.EQ.1) GOTO 31
      TA=C
      DO 30 K=2,J
      JK=J-K+1
      B1=D(J,JK)
      W=C-V
      U=W/(B1-ONE)
      C=B1*U
      V=DT(I,K)
      DT(I,K)=U
30    TA=U+TA
      US(I)=TA
      TA=U/S(I)
      ERR=ERR+TA*TA
31    CONTINUE
      IF(J.EQ.1) GOTO 35
C ERROR (SCALED ROOT MEAN SQUARE)
      ERR=DSQRT(ERR/FN)
CSP       ERR=SQRT(ERR/FN)
      KONV=0
      IF(ERR.LT.EPS) KONV=1
      ERR=ERR/EPH
C
C ORDER CONTROL
      K=J-1
      L=J+K
      FC=ERR**(ONE/DFLOAT(L))
CSP       FC=ERR**(ONE/FLOAT(L))
      IF(FC.LT.FCM) FC=FCM
C  OPTIMAL ORDER DETERMINATION
      OMJ=FC*A(J)
      IF(J.GT.2.AND.OMJ*ONE1.GT.OMJO.OR.K.GT.JOH) GOTO 32
      KO=K
      JO=J
      OMJO=OMJ
      FCO=FC
32    IF(J.LT.KOH.AND.NSTEP.GT.0) GOTO 35
      IF(KONV.EQ.0) GOTO 33
      IF(KO.LT.K.OR.INCR(J).LT.0) GOTO 40
C  POSSIBLE INCREASE OF ORDER
      IF(NRED(KO).GT.0) NRED(KO)=NRED(KO)-1
      IF(J.GE.JM) GOTO 40
      FC=FCO/AL(J,K)
      IF(FC.LT.FCM) FC=FCM
      J1=J+1
      IF(A(J1)*FC*ONE1.GT.OMJO) GOTO 40
      FCO=FC
      KO=JO
      JO=JO+1
      GOTO 40
C
C
C  CONVERGENCE MONITOR
33    RED=ONE/FCO
      JK=KM
      IF(JOH.LT.KM) JK=JOH
      IF(K.GE.JK) GOTO 36
      IF(KO.LT.KOH) RED=AL(KOH,KO)/FCO
34    IF(AL(JK,KO).LT.FCO) GOTO 36
35    CONTINUE
C
C STEPSIZE REDUCTION
36    RED=RED*SAFE
      H=H*RED
      IF(NSTEP.EQ.0) GOTO 38
      NRED(KOH)=NRED(KOH)+1
      DO 37 L=KOH,KM
37    INCR(L)=-2-NRED(KOH)
38    JRED=JRED+1
      IF(KFLAG.GE.1) WRITE(LOUT,1002) JRED,RED,KOH
      IF(JRED.GT.JRMAX) GOTO 52
      GOTO 20
C
C  PREPARATIONS FOR NEXT BASIC INTEGRATION STEP
40    T=TN
      H1=TEND-T
      DO  41 I=1,NF
      Y(I)=UH(I)+US(I)
41    Y(NF+I)=US(NF+I)+ZH(I)
      NSTEP=NSTEP+1
      IF(NSTEP.GT.NSTMAX) GO TO 51
C
C STEPSIZE PREDICTION
      IF(FCO.NE.FCM) HR=H
      H=H/FCO
      KOH=KO
      JOH=KOH+1
      IF(DABS(H).LE.DABS(T)*EPMACH) GO TO 53
CSP       IF(ABS(H).LE.ABS(T)*EPMACH) GO TO 53
      HR=H
42    CONTINUE
      IF(DABS(H1).LE.TEPS) GOTO 50
CSP       IF(ABS(H1).LE.TEPS) GOTO 50
      Q=H1/H
      IF(Q.GT.EPMACH) GOTO 15
C
C  SOLUTION EXIT
50    H=HR
      IF(KFLAG.GT.0) WRITE(LOUT,1009) NSTEP,NFCN,T,K,KOH
      IF(KFLAG.GT.1) WRITE(LOUT,1000) NSTEP,NFCN,T,(Y(I),I=1,N)
      HMAX=HMAXU
      RETURN
C
C  FAIL EXIT
51    IF(KFLAG.GE.1)WRITE(LOUT,1008) NSTMAX
      KFLAG=-1
      GOTO 55
52    IF(KFLAG.GE.1) WRITE(LOUT,1010)JRMAX
      KFLAG=-2
      GOTO 55
53    IF(KFLAG.GE.1) WRITE(LOUT,1004)
      KFLAG=-3
      GOTO 55
54    IF(KFLAG.GE.1)WRITE(LOUT,1011)
      KFLAG=-4
55    H=ZERO
      HMAX=HMAXU
      IF(KFLAG.GE.1) WRITE(LOUT,1009) NSTMAX,NFCN,T,K,KO
      IF(KFLAG.GT.1) WRITE(LOUT,1000) NSTMAX,NFCN,T,(Y(I),I=1,N)
      RETURN
C
1000  FORMAT(1H ,2I9,5D20.11,/,(1H ,38X,4D20.11))
CSP 1000  FORMAT(1H ,2I9,5E20.6,/,(1H ,38X,4E20.6))
1001  FORMAT(1H0,27H  KEPLEX     REL.PREC. EPS ,D10.3,8HMAX.COL.,I3)
CSP 1001  FORMAT(1H0,27H  KEPLEX     REL.PREC. EPS ,E10.3,8HMAX.COL.,I3)
1002  FORMAT(1H ,I3,17HREDUCTION FACTOR ,D10.3,I9,/)
CSP 1002  FORMAT(1H ,I3,17HREDUCTION FACTOR ,E10.3,I9,/)
1004  FORMAT(/,40H0  STEPSIZE REDUCTION FAILED TO SUCCEED ,//)
1005  FORMAT(//,5X,4HSTEP,3X,7HF-CALLS,8X,1HX,19X,7HY1(X)..,//)
1008  FORMAT(18H0MORE THAN NSTMAX=,I6,18H INTEGRATION STEPS,//)
1009  FORMAT(1H ,2I9,D20.11,I9,I6,/)
CSP 1009  FORMAT(1H ,2I9,E20.6,I9,I6,/)
1010  FORMAT(17H0MORE THAN JRMAX=,I3,29H STEPSIZE REDUCTIONS PER STEP,/)
1011  FORMAT(29H0 N ODD,KEPLEX NOT APPLICABLE,//)
C
C
C
C  END KEPLEX
C
      END
