      SUBROUTINE DIFEX2 (N,FCN,T,Y,TEND,EPS,HMAX,H,KFLAG)
C
C* Begin Prologue DIFEX2
C
C  ---------------------------------------------------------------------
C
C* Title
C
C    Explicit extrapolation integrator for non-stiff systems of
C    second-order differential equations with the first derivative
C    absent in the right-hand side (based on Stoermer discretization).
C
C    Special version adapted for possible call in the multiple shooting
C    code BVPSOL.
C
C* Written by        P. Deuflhard, U. Nowak, U. Poehle
C* Purpose           Solution of systems of initial value problems
C* Method            Explicit mid-point rule discretization with
C                    h**2-extrapolation
C* Category          i1a1c2. - System of nonstiff second order
C                              differential equations
C* Keywords          extrapolation, ODE, explicit mid-point rule,
C                    non-stiff
C* Version           1.0 , August 1982
C* Latest Change     Februar 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  REFERENCES:
C
C /1/ W.B. GRAGG:
C     ON EXTRAPOLATION ALGORITHMS FOR ORDINARY INITIAL VALUE PROBLEMS
C     SIAM J. NUMER. ANAL. 2, 384-404 (1965)
C
C /2/ P.DEUFLHARD:
C     ORDER AND STEPSIZE CONTROL IN EXTRAPOLATION METHODS
C     NUMER. MATH. 41, 373-398 (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 VALUES INCLUDING DERIVATIVES
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 DIFEX2 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),YM(12),DY2( 6),DZ2( 6),S(12)
CSP       REAL         Y(N),YM(12),DY2( 6),DZ2( 6),S(12)
C                      Y(N),YM( N),DY2(NF),DZ2(NF),S( N)
      DOUBLE PRECISION DT(12,13),D(13,13),A(13),AL(13,13)
CSP       REAL         DT(12,13),D(13,13),A(13),AL(13,13)
C                      DT( N,JM),D(JM,JM),A(JM),AL(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 FN,G,H,HALF,HMAX,HMAXU,HR,H1,OMJ,OMJO,ONE
CSP       REAL FN,G,H,HALF,HMAX,HMAXU,HR,H1,OMJ,OMJO,ONE
      DOUBLE PRECISION RED,RO,SAFE,SMALL,TA,TEN,U,U1,V,W,T,TEND
CSP       REAL RED,RO,SAFE,SMALL,TA,TEN,U,U1,V,W,T,TEND
      DOUBLE PRECISION B2,FMIN,ONE1,Q,TEPS,TN,ZERO
CSP       REAL B2,FMIN,ONE1,Q,TEPS,TN,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.50D0/
CSP       DATA ZERO/0.E0/,FMIN/1.E-2/,RO/0.25E0/,HALF/0.50E0/
      DATA ONE/1.D0/,ONE1/1.01D0/,TEN/1.D1/,SAFE/0.5D0/
CSP       DATA ONE/1.E0/,ONE1/1.01E0/,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/
C
C  SET MACHINE CONSTANTS
      CALL ZIBCONST(EPMACH, SMALL)
C
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
      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.GT.0) WRITE(LOUT,1001) EPS,KM
      IF(KFLAG.GT.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
C
C  SCALING
      DO 18 I=1,N
      U=DABS(Y(I))
CSP       U=ABS(Y(I))
      IF(U.LT.EPMACH) U=ONE
18    S(I)=U
      CALL FCN(NF,T,Y,DZ2)
      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  STOERMER DISCRETIZATION
      DO 35 J=1,JM
      M=NJ(J)
      B=H/DFLOAT(M)
CSP       B=H/FLOAT(M)
      G=B*HALF
      B2=B*B
      DO 21 I=1,NF
      NFI=NF+I
      V=Y(NFI)+G*DZ2(I)
      V=V*B
      YM(NFI)=V
21    YM(I)=Y(I)+V
      M=M-1
      IF(M.EQ.0) GOTO 23
      DO 22 K=1,M
      CALL FCN(NF,T+DFLOAT(K)*B,YM,DY2)
CSP       CALL FCN(NF,T+FLOAT(K)*B,YM,DY2)
      NFCN=NFCN+1
      DO 22 I=1,NF
      NFI=NF+I
      V=YM(NFI)+B2*DY2(I)
      YM(NFI)=V
      YM(I)=YM(I)+V
22    CONTINUE
23    CALL FCN(NF,TN,YM,DY2)
      NFCN=NFCN+1
      DO 24 I=1,NF
      NFI=NF+I
24    YM(NFI)=YM(NFI)/B+G*DY2(I)
C
C  EXTRAPOLATION
C
      ERR=ZERO
      DO 31 I=1,N
      V=DT(I,1)
      C=YM(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
      YM(I)=TA
C PREPARATIONS FOR RE-SCALING
      TA=DABS(TA)
CSP       TA=ABS(TA)
      IF(TA.LT.S(I)) TA=S(I)
      U=U/TA
      ERR=ERR+U*U
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.GT.0) 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,N
41    Y(I)=YM(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.1) 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.GT.0)WRITE(LOUT,1008) NSTMAX
      KFLAG=-1
      GOTO 55
52    IF(KFLAG.GT.0) WRITE(LOUT,1010)JRMAX
      KFLAG=-2
      GOTO 55
53    IF(KFLAG.GT.0) WRITE(LOUT,1004)
      KFLAG=-3
      GOTO 55
54    IF(KFLAG.GT.0)WRITE(LOUT,1011)
      KFLAG=-4
55    H=ZERO
      HMAX=HMAXU
      IF(KFLAG.GT.0) 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  DIFEX2     REL.PREC. EPS ,D10.3,8HMAX.COL.,I3)
CSP 1001  FORMAT(1H0,27H  DIFEX2     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,DIFEX2 NOT APPLICABLE,//)
C
C
C
C  END DIFEX2
C
      END
