SUBROUTINE ODEX2(N,FCN,X,Y,YP,XEND,H,
     &                 RTOL,ATOL,ITOL,
     &                 SOLOUT,IOUT,
     &                 WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
C ----------------------------------------------------------
C     NUMERICAL SOLUTION OF A SYSTEM OF SECOND 0RDER
C     ORDINARY DIFFERENTIAL EQUATIONS  Y''=F(X,Y).
C     THIS IS AN EXTRAPOLATION-ALGORITHM, BASED ON
C     THE STOERMER RULE (WITH STEPSIZE CONTROL
C     ORDER SELECTION AND DENSE OUTPUT).
C     
C     AUTHORS: E. HAIRER AND G. WANNER
C              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
C              CH-1211 GENEVE 24, SWITZERLAND 
C              E-MAIL:  Ernst.Hairer@math.unige.ch
C                       Gerhard.Wanner@math.unige.ch
C     
C     THIS CODE IS DESCRIBED IN SECTION II.14 OF THE BOOK:
C         E. HAIRER, S.P. NORSETT AND G. WANNER, SOLVING ORDINARY
C         DIFFERENTIAL EQUATIONS I. NONSTIFF PROBLEMS. 2ND EDITION.
C         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
C         SPRINGER-VERLAG (1993) 
C
C     VERSION SEPTEMBER 30, 1995
C
C     INPUT PARAMETERS  
C     ----------------  
C     N           DIMENSION OF THE SYSTEM 
C
C     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
C                 VALUE OF F(X,Y):
C                    SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR)
C                    DOUBLE PRECISION X,Y(N),F(N)
C                    F(1)=...   ETC.
C
C     X           INITIAL X-VALUE
C
C     Y(N)        INITIAL VALUES FOR Y
C
C     YP(N)       INITIAL VALUES FOR Y'
C
C     XEND        FINAL X-VALUE
C
C     H           INITIAL STEP SIZE GUESS;
C                 USUALLY 1.D-1 OR 1.D-3, IS GOOD. THIS CHOICE IS NOT 
C                 VERY IMPORTANT, THE CODE QUICKLY ADAPTS ITS STEPSIZE.
C                 WHEN YOU ARE NOT SURE, THEN STUDY THE CHOSEN VALUES
C                 FOR A FEW STEPS IN SUBROUTINE "SOLOUT".
C                 (IF H=0.D0, THE CODE PUTS H=1.D-4).
C
C     RTOL,ATOL   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY CAN
C                 BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH 2*N.
C
C     ITOL        SWITCH FOR RTOL AND ATOL:
C                   ITOL=0: BOTH RTOL AND ATOL ARE SCALARS.
C                     THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF
C                     Y(I)  BELOW  RTOL*ABS(Y(I))+ATOL
C                     YP(I) BELOW  RTOL*ABS(YP(I))+ATOL
C                   ITOL=1: BOTH RTOL AND ATOL ARE VECTORS.
C                     THE CODE KEEPS THE LOCAL ERROR OF
C                     Y(I)  BELOW  RTOL(I)*ABS(Y(I))+ATOL(I).
C                     YP(I) BELOW  RTOL(I+N)*ABS(YP(I))+ATOL(I+N).
C
C     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
C                 NUMERICAL SOLUTION DURING INTEGRATION. 
C                 IF IOUT>=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
C                 SUPPLY A DUMMY SUBROUTINE IF IOUT=0. 
C                 IT MUST HAVE THE FORM
C                    SUBROUTINE SOLOUT (NR,XOLD,X,Y,YP,N,CON,NCON,ICOMP,ND,
C                                       RPAR,IPAR,IRTRN)
C                    DIMENSION Y(N),YP(N),CON(NCON),ICOMP(ND)
C                    ....  
C                 SOLOUT FURNISHES THE SOLUTIONS "Y, YP" AT THE NR-TH
C                    GRID-POINT "X" (THEREBY THE INITIAL VALUE IS
C                    THE FIRST GRID-POINT).
C                 "XOLD" IS THE PRECEEDING GRID-POINT.
C                 "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
C                    IS SET <0, ODEX2 WILL RETURN TO THE CALLING PROGRAM.
C           
C          -----  CONTINUOUS OUTPUT (IF IOUT=2): -----
C                 DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
C                 FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH
C                 THE DOUBLE PRECISION FUNCTION
C                    >>>   CONTX2(I,S,CON,NCON,ICOMP,ND)   <<<
C                 WHICH PROVIDES AN APPROXIMATION TO THE I-TH
C                 COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
C                 S SHOULD LIE IN THE INTERVAL [XOLD,X].
C           
C     IOUT        SWITCH FOR CALLING THE SUBROUTINE SOLOUT:
C                    IOUT=0: SUBROUTINE IS NEVER CALLED
C                    IOUT=1: SUBROUTINE IS USED FOR OUTPUT
C                    IOUT=2: DENSE OUTPUT IS PERFORMED IN SOLOUT
C
C     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
C                 SERVES AS WORKING SPACE FOR ALL VECTORS.
C                 "LWORK" MUST BE AT LEAST
C                    N*(2*KM+6)+5*KM+20+(KM*(2*KM+5)+6)*NRDENS
C                 WHERE NRDENS=IWORK(8) (SEE BELOW) AND
C                        KM=9                IF IWORK(2)=0
C                        KM=IWORK(2)         IF IWORK(2).GT.0
C                 WORK(1),...,WORK(20) SERVE AS PARAMETERS
C                 FOR THE CODE. FOR STANDARD USE, SET THESE
C                 PARAMETERS TO ZERO BEFORE CALLING.
C
C     LWORK       DECLARED LENGTH OF ARRAY "WORK".
C
C     IWORK       INTEGER WORKING SPACE OF LENGTH "LIWORK".
C                 "LIWORK" MUST BE AT LEAST
C                               KM+20+NRDENS 
C                 IWORK(1),...,IWORK(20) SERVE AS PARAMETERS
C                 FOR THE CODE. FOR STANDARD USE, SET THESE
C                 PARAMETERS TO ZERO BEFORE CALLING.
C
C     LIWORK      DECLARED LENGTH OF ARRAY "IWORK".
C
C     RPAR, IPAR  REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH  
C                 CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING
C                 PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES. 
C
C-----------------------------------------------------------------------
C 
C     SOPHISTICATED SETTING OF PARAMETERS
C     -----------------------------------
C              SEVERAL PARAMETERS (WORK(1),...,IWORK(1),...) ALLOW
C              TO ADAPT THE CODE TO THE PROBLEM AND TO THE NEEDS OF
C              THE USER. FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES.
C
C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 2.3D-16.
C
C    WORK(2)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
C
C    WORK(3)   STEP SIZE IS REDUCED BY FACTOR WORK(3), IF DURING THE
C              COMPUTATION OF THE EXTRAPOLATION TABLEAU DIVERGENCE
C              IS OBSERVED; DEFAULT 0.5.
C
C    WORK(4), WORK(5)   PARAMETERS FOR STEP SIZE SELECTION
C              THE NEW STEP SIZE FOR THE J-TH DIAGONAL ENTRY IS
C              CHOSEN SUBJECT TO THE RESTRICTION
C                 FACMIN/WORK(5) <= HNEW(J)/HOLD <= 1/FACMIN
C              WHERE FACMIN=WORK(4)**(1/(2*J-1)) 
C              DEFAULT VALUES: WORK(4)=0.02D0, WORK(5)=4.D0
C
C    WORK(6), WORK(7)   PARAMETERS FOR THE ORDER SELECTION
C              STEP SIZE IS DECREASED IF    W(K-1) <= W(K)*WORK(6)
C              STEP SIZE IS INCREASED IF    W(K) <= W(K-1)*WORK(7)
C              DEFAULT VALUES: WORK(6)=0.8D0, WORK(7)=0.9D0
C
C    WORK(8), WORK(9)   SAFETY FACTORS FOR STEP CONTROL ALGORITHM
C             HNEW=H*WORK(9)*(WORK(8)*TOL/ERR)**(1/(J-1))
C             DEFAULT VALUES: WORK(8)=0.65D0,
C                        WORK(9)=0.94D0  IF "HOPE FOR CONVERGENCE"
C
C    IWORK(1)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
C              THE DEFAULT VALUE (FOR IWORK(1)=0) IS 10000.
C
C    IWORK(2)  THE MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION 
C              TABLE. THE DEFAULT VALUE (FOR IWORK(2)=0) IS 9.
C              IF IWORK(2).NE.0 THEN IWORK(2) SHOULD BE .GE.3.
C
C    IWORK(3)  SWITCH FOR THE STEP SIZE SEQUENCE (EVEN NUMBERS ONLY)
C              IF IWORK(3).EQ.1 THEN 2,4,6,8,10,12,14,16,...
C              IF IWORK(3).EQ.2 THEN 2,4,8,12,16,20,24,28,...
C              IF IWORK(3).EQ.3 THEN 2,4,6,8,12,16,24,32,...
C              IF IWORK(3).EQ.4 THEN 2,6,10,14,18,22,26,30,...
C              THE DEFAULT VALUE IS IWORK(3)=1 IF IOUT.LE.1;
C              THE DEFAULT VALUE IS IWORK(3)=4 IF IOUT.GE.2.
C 
C    IWORK(6)  IF  IWORK(6)=0  ERROR ESTIMATOR IN THE DENSE
C              OUTPUT FORMULA IS ACTIVATED. IT CAN BE SUPPRESSED
C              BY PUTTING IWORK(6)=1.
C              DEFAULT IWORK(6)=0  (IF IOUT.GE.2).
C
C    IWORK(7)  DETERMINES THE DEGREE OF INTERPOLATION FORMULA 
C              MU = 2 * KAPPA - IWORK(7) + 1
C              IWORK(7) SHOULD LIE BETWEEN 1 AND 8
C              DEFAULT IWORK(7)=6  (IF IWORK(7)=0).
C
C    IWORK(8)  = NRDENS = NUMBER OF COMPONENTS, FOR WHICH DENSE OUTPUT
C              IS REQUIRED
C
C    IWORK(21),...,IWORK(NRDENS+20) INDICATE THE COMPONENTS, FOR WHICH
C              DENSE OUTPUT IS REQUIRED. THEY NEED NOT BE SET IF
C              NRDENS = N .
C
C----------------------------------------------------------------------C
C     OUTPUT PARAMETERS 
C     ----------------- 
C     X           X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED
C                 (AFTER SUCCESSFUL RETURN X=XEND).
C
C     Y(N)        NUMERICAL SOLUTION AT X
C 
C     YP(N)       NUMERICAL DERIVATIVE OF SOLUTION AT X
C 
C     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
C
C     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
C                   IDID=1  COMPUTATION SUCCESSFUL,
C                   IDID=-1 COMPUTATION UNSUCCESSFUL.
C
C   IWORK(17)  NFCN    NUMBER OF FUNCTION EVALUATIONS
C   IWORK(18)  NSTEP   NUMBER OF COMPUTED STEPS
C   IWORK(19)  NACCPT  NUMBER OF ACCEPTED STEPS
C   IWORK(20)  NREJCT  NUMBER OF REJECTED STEPS (DUE TO ERROR TEST),
C                      (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED)
C-----------------------------------------------------------------------
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C          DECLARATIONS 
C *** *** *** *** *** *** *** *** *** *** *** *** ***
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(N),YP(N),WORK(LWORK),IWORK(LIWORK)
      DIMENSION ATOL(*),RTOL(*),RPAR(*),IPAR(*)
      LOGICAL ARRET
      EXTERNAL FCN,SOLOUT
C *** *** *** *** *** *** ***
C        SETTING THE PARAMETERS 
C *** *** *** *** *** *** *** 
      NFCN=0
      NSTEP=0
      NACCPT=0
      NREJCT=0
      ARRET=.FALSE.
C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
      IF(IWORK(1).EQ.0)THEN
         NMAX=10000
      ELSE
         NMAX=IWORK(1)
         IF(NMAX.LE.0)THEN
            WRITE(6,*)' WRONG INPUT IWORK(1)=',IWORK(1)
            ARRET=.TRUE.
         END IF
      END IF
C -------- KM     MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION 
      IF(IWORK(2).EQ.0)THEN
         KM=9
      ELSE
         KM=IWORK(2)
         IF(KM.LE.2)THEN
            WRITE(6,*)' CURIOUS INPUT IWORK(2)=',IWORK(2)
            ARRET=.TRUE.
         END IF
      END IF
C -------- NSEQU     CHOICE OF STEP SIZE SEQUENCE
      NSEQU=IWORK(3)
      IF(IWORK(3).EQ.0.AND.IOUT.LE.1) NSEQU=1
      IF(IWORK(3).EQ.0.AND.IOUT.GE.2) NSEQU=4
      IF(NSEQU.LE.0.OR.NSEQU.GE.5)THEN
         WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
         ARRET=.TRUE.
      END IF 
      IF (NSEQU.LE.3.AND.IOUT.GE.2) THEN
         WRITE(6,*)' IWORK(3) NOT COMPATIBLE WITH IOUT'
         ARRET=.TRUE.
      END IF 
C -------- IDERR  PARAMETER FOR ERROR ESTIMATION IN DENSE OUTPUT
      IF(IWORK(6).EQ.0)THEN
         IF(IOUT.LE.1) IDERR=1
         IF(IOUT.GE.2) IDERR=0
      ELSE
         IDERR=IWORK(6)
         IF(IOUT.LE.1)THEN
            WRITE(6,*)' ERROR ESTIMATION IN DENSE OUTPUT',
     &                ' NOT POSSIBLE, WRONG IWORK(6)=',IWORK(6)
            ARRET=.TRUE.
         END IF
      END IF
C -------- MUDIF
      IF(IWORK(7).EQ.0)THEN
         MUDIF=6
      ELSE
         MUDIF=IWORK(7)
         IF(MUDIF.LE.0.OR.MUDIF.GE.9)THEN
            WRITE(6,*)' WRONG INPUT IWORK(7)=',IWORK(7)
            ARRET=.TRUE.
         END IF
      END IF
C -------- NRDENS   NUMBER OF DENSE OUTPUT COMPONENTS
      NRDENS=IWORK(8)
      IF(NRDENS.LT.0.OR.NRDENS.GT.N)THEN
         WRITE(6,*)' CURIOUS INPUT IWORK(8)=',IWORK(8)
         ARRET=.TRUE.
      END IF
      IF (NRDENS.EQ.N) THEN
           DO 17 I=1,NRDENS
  17       IWORK(20+I)=I 
      END IF
C -------- UROUND   SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0  
      IF(WORK(1).EQ.0.D0)THEN
         UROUND=2.3D-16
      ELSE
         UROUND=WORK(1)
         IF(UROUND.LE.1.D-35.OR.UROUND.GE.1.D0)THEN
            WRITE(6,*)' WHICH MACHINE DO YOU HAVE? YOUR UROUND WAS:'
     &                            ,WORK(1)
            ARRET=.TRUE.
         END IF
      END IF
C -------- MAXIMAL STEP SIZE
      IF(WORK(2).EQ.0.D0)THEN
         HMAX=XEND-X
      ELSE
         HMAX=ABS(WORK(2))
      END IF
C -------- STEP SIZE REDUCTION FACTOR
      IF(WORK(3).EQ.0.D0)THEN
         SAFE3=0.5D0
      ELSE
         SAFE3=WORK(3)
         IF(SAFE3.LE.UROUND.OR.SAFE3.GE.1.D0)THEN
            WRITE(6,*)' CURIOUS INPUT WORK(3)=',WORK(3)
            ARRET=.TRUE.
         END IF
      END IF
C -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
      IF(WORK(4).EQ.0.D0)THEN
         FAC1=0.02D0
      ELSE
         FAC1=WORK(4)
      END IF
      IF(WORK(5).EQ.0.D0)THEN
         FAC2=4.0D0
      ELSE
         FAC2=WORK(5)
      END IF
C -------  FAC3, FAC4   PARAMETERS FOR THE ORDER SELECTION
      IF(WORK(6).EQ.0.D0)THEN
         FAC3=0.8D0
      ELSE
         FAC3=WORK(6)
      END IF
      IF(WORK(7).EQ.0.D0)THEN
         FAC4=0.9D0
      ELSE
         FAC4=WORK(7)
      END IF
C ------- SAFE1, SAFE2 SAFETY FACTORS FOR STEP SIZE PREDICTION
      IF(WORK(8).EQ.0.D0)THEN
         SAFE1=0.65D0
      ELSE
         SAFE1=WORK(8)
      END IF
      IF(WORK(9).EQ.0.D0)THEN
         SAFE2=0.94D0
      ELSE
         SAFE2=WORK(9)
      END IF
C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
      LFSAFE=2*KM*KM+KM 
      IEDY=21
      IEYH1=IEDY+N
      IEYH2=IEYH1+N
      IEDZ=IEYH2+N
      IESCAL=IEDZ+N
      IESCP=IESCAL+N
      IET=IESCP+N
      IETP=IET+KM*N
      IEFS=IETP+KM*N
      IEYS=IEFS+LFSAFE*NRDENS
      IEYSP=IEYS+KM*NRDENS
      IEHH=IEYSP+KM*NRDENS
      IEW=IEHH+KM
      IEA=IEW+KM 
      IEFAC=IEA+KM
C ------ TOTAL STORAGE REQUIREMENT -----------
      IECO=IEFAC+2*KM
      ISTORE=IECO+(2*KM+6)*NRDENS-1
      IF(ISTORE.GT.LWORK)THEN
         WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
         ARRET=.TRUE.
      END IF
C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
      ICOM=21
      IENJ=ICOM+NRDENS
C --------- TOTAL REQUIREMENT ---------------
      ISTORE=IENJ+KM-1
      IF(ISTORE.GT.LIWORK)THEN
         WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
         ARRET=.TRUE.
      END IF
C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
      IF (ARRET) THEN
         IDID=-1
         RETURN
      END IF
C -------- CALL TO CORE INTEGRATOR ------------
      NRD=MAX(1,NRDENS)  
      NCOM=MAX(1,(2*KM+6)*NRDENS)
      CALL ODX2CO(N,FCN,X,Y,YP,XEND,HMAX,H,RTOL,ATOL,ITOL,KM,
     &   SOLOUT,IOUT,IDID,NMAX,UROUND,WORK(IEDY),WORK(IEYH1),
     &   WORK(IEYH2),WORK(IEDZ),WORK(IESCAL),WORK(IESCP),WORK(IEFS),
     &   WORK(IEYS),WORK(IEYSP),WORK(IET),WORK(IETP),WORK(IEHH),
     &   WORK(IEW),WORK(IEA),WORK(IECO),NCOM,IWORK(ICOM),IWORK(IENJ),
     &   NSEQU,LFSAFE,SAFE1,SAFE2,SAFE3,FAC1,FAC2,FAC3,FAC4,IDERR,
     &   WORK(IEFAC),MUDIF,NRD,RPAR,IPAR,NFCN,NSTEP,NACCPT,NREJCT)
      IWORK(17)=NFCN
      IWORK(18)=NSTEP
      IWORK(19)=NACCPT
      IWORK(20)=NREJCT
C ----------- RETURN -----------
      RETURN
      END
C
C
C
C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
C
      SUBROUTINE ODX2CO(N,FCN,X,Y,YP,XEND,HMAX,H,RTOL,ATOL,ITOL,KM,
     &   SOLOUT,IOUT,IDID,NMAX,UROUND,DY,YH1,YH2,DZ,SCAL,SCALP,FSAFE,
     &   YSAFE,YSAFEP,T,TP,HH,W,A,DENS,NCOM,ICOMP,NJ,NSEQU,LFSAFE,
     &   SAFE1,SAFE2,SAFE3,FAC1,FAC2,FAC3,FAC4,IDERR,ERRFAC,MUDIF,
     &   NRD,RPAR,IPAR,NFCN,NSTEP,NACCPT,NREJCT)
C ----------------------------------------------------------
C     CORE INTEGRATOR FOR ODEX2
C     PARAMETERS SAME AS IN ODEX2 WITH WORKSPACE ADDED 
C ---------------------------------------------------------- 
C         DECLARATIONS 
C ---------------------------------------------------------- 
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       LOGICAL REJECT,LAST,ATOV
       EXTERNAL FCN
       DIMENSION Y(N),YP(N),DY(N),YH1(N),YH2(N),DZ(N),SCAL(N),SCALP(N) 
       DIMENSION T(KM,N),TP(KM,N),NJ(KM),HH(KM),W(KM),A(KM)
       DIMENSION RTOL(*),ATOL(*),RPAR(*),IPAR(*)
       DIMENSION FSAFE(LFSAFE,NRD),YSAFE(KM,NRD),YSAFEP(KM,NRD)   
       DIMENSION ERRFAC(2*KM),DENS(NCOM),ICOMP(NRD)
       COMMON /CONOD2/XOLDD,HHH,KMIT
C --- DEFINE THE STEP SIZE SEQUENCE
       IF (NSEQU.EQ.1) THEN
           DO 1 I=1,KM
   1       NJ(I)=2*I
       END IF
       IF (NSEQU.EQ.2) THEN 
           NJ(1)=2
           DO 2 I=2,KM
   2       NJ(I)=4*I-4
       END IF
       IF (NSEQU.EQ.3) THEN 
           NJ(1)=2
           NJ(2)=4
           NJ(3)=6
           DO 11 I=4,KM
   11      NJ(I)=2*NJ(I-2)
       END IF
       IF (NSEQU.EQ.4) THEN
           DO 3 I=1,KM
   3       NJ(I)=4*I-2
       END IF
C --- DEFINE THE A(I) FOR ORDER SELECTION
       A(1)=1.D0+NJ(1)/2
       DO 4 I=2,KM
   4   A(I)=A(I-1)+NJ(I)/2
C --- INITIAL SCALING
       IF (ITOL.EQ.0) THEN
         DO 6 I=1,N
         SCAL(I)=ATOL(1)+RTOL(1)*ABS(Y(I))
  6      SCALP(I)=ATOL(1)+RTOL(1)*ABS(YP(I))
       ELSE
       DO 8 I=1,N
         SCAL(I)=ATOL(I)+RTOL(I)*ABS(Y(I))
  8      SCALP(I)=ATOL(I+N)+RTOL(I+N)*ABS(YP(I))
       END IF
C --- INITIAL PREPARATIONS
       K=MAX(3,MIN(KM-1,INT(-LOG10(RTOL(1))*0.6D0+1.5D0)))
       IF (H.LE.UROUND) H=1.D-4 
       H=MIN(H,HMAX,XEND-X)
       IF (IOUT.GE.1) THEN
         IF (IOUT.GE.2) THEN
            DO 9 MU=1,KM*2
            ERRX=SQRT(MU/(MU+6.D0))*0.5D0
            PROD=(1.5D0/(MU+6.D0))**3
            DO 7 J=1,MU
   7        PROD=PROD*ERRX/J
   9        ERRFAC(MU)=PROD 
         END IF
         NNRD=NRD
         IRTRN=0
         XOLD=X
         CALL SOLOUT (NACCPT+1,XOLD,X,Y,YP,N,DENS,NCON,ICOMP,NRD,
     &                 RPAR,IPAR,IRTRN)
         IF (IRTRN.LT.0) GOTO 120
       END IF 
       ERR=0.D0
       ERROLD=1.D10  
       HOPTDE=HMAX
       W(1)=0.D0  
       REJECT=.FALSE.
       LAST=.FALSE.
  10   ATOV=.FALSE.
C --- IS XEND REACHED IN THE NEXT STEP?
       H=MIN(H,HMAX,HOPTDE)
       IF (X+H.GE.XEND) THEN
          LAST=.TRUE.
          H=XEND-X
       END IF
       IF (NSTEP.EQ.0.OR.IOUT.NE.2) CALL FCN(N,X,Y,DZ,RPAR,IPAR)
       NFCN=NFCN+1
C --- THE FIRST AND LAST STEP 
       IF (NSTEP.EQ.0.OR.LAST) THEN 
          IPT=0
          NSTEP=NSTEP+1 
          DO 20 J=1,K
          KC=J 
          CALL STOERM(J,X,Y,YP,H,HMAX,N,FCN,DY,YH1,YH2,DZ,T,TP,NJ,HH,W,
     1      ERR,FAC,A,SAFE1,UROUND,FAC1,FAC2,SAFE2,SCAL,SCALP,ATOV,
     2      SAFE3,REJECT,KM,RTOL,ATOL,ITOL,ERROLD,FSAFE,
     3      LFSAFE,IOUT,IPT,YSAFE,YSAFEP,ICOMP,NRD,RPAR,IPAR,NFCN)
          IF (ATOV) GO TO 10
  20      IF (J.GT.1.AND.ERR.LE.1.D0) GO TO 60
          GO TO 55
       END IF
C --- BASIC INTEGRATION STEP  
  30   CONTINUE
       IPT=0
       NSTEP=NSTEP+1
       IF (NSTEP.GE.NMAX) GO TO 120 
       KC=K-1
       DO 40 J=1,KC
       CALL STOERM(J,X,Y,YP,H,HMAX,N,FCN,DY,YH1,YH2,DZ,T,TP,NJ,HH,W,
     1   ERR,FAC,A,SAFE1,UROUND,FAC1,FAC2,SAFE2,SCAL,SCALP,ATOV,SAFE3,
     2   REJECT,KM,RTOL,ATOL,ITOL,ERROLD,FSAFE,LFSAFE,
     3      IOUT,IPT,YSAFE,YSAFEP,ICOMP,NRD,RPAR,IPAR,NFCN)
       IF (ATOV) GO TO 10
  40   CONTINUE
C --- CONVERGENCE MONITOR
       IF (K.EQ.2.OR.REJECT) GO TO 50
       IF (ERR.LE.1.D0) GO TO 60
       IF (ERR.GT.((NJ(K+1)*NJ(K))/4.D0)**2) GO TO 100  
  50   CONTINUE
       CALL STOERM(K,X,Y,YP,H,HMAX,N,FCN,DY,YH1,YH2,DZ,T,TP,NJ,HH,W,
     1   ERR,FAC,A,SAFE1,UROUND,FAC1,FAC2,SAFE2,SCAL,SCALP,ATOV,SAFE3,
     2   REJECT,KM,RTOL,ATOL,ITOL,ERROLD,FSAFE,LFSAFE,
     3      IOUT,IPT,YSAFE,YSAFEP,ICOMP,NRD,RPAR,IPAR,NFCN)
       IF (ATOV) GO TO 10
       KC=K 
       IF (ERR.LE.1.D0) GO TO 60
C --- HOPE FOR CONVERGENCE IN LINE K+1
  55   CONTINUE
       IF (ERR.GT.(NJ(K+1)/2.D0)**2) GO TO 100  
       KC=K+1
       CALL STOERM(KC,X,Y,YP,H,HMAX,N,FCN,DY,YH1,YH2,DZ,T,TP,NJ,HH,W,
     1   ERR,FAC,A,SAFE1,UROUND,FAC1,FAC2,SAFE2,SCAL,SCALP,ATOV,SAFE3,
     2   REJECT,KM,RTOL,ATOL,ITOL,ERROLD,FSAFE,LFSAFE,
     3      IOUT,IPT,YSAFE,YSAFEP,ICOMP,NRD,RPAR,IPAR,NFCN)
       IF (ATOV) GO TO 10
       IF (ERR.GT.1.D0) GO TO 100
C --- STEP IS ACCEPTED  
  60   XOLD=X
       X=X+H
       IF (IOUT.GE.2) THEN
C ---  KMIT = MU OF THE PAPER
           KMIT=MAX(1,2*KC-MUDIF+1)
           DO 478 I=1,N
 478       YH1(I)=T(1,I)
           CALL FCN(N,X,YH1,YH2,RPAR,IPAR)
           H2=H*H
           DO 69 I=1,NRD 
           DENS(I)=Y(ICOMP(I))  
           DENS(NRD+I)=H*YP(ICOMP(I))  
           DENS(2*NRD+I)=H2*DZ(ICOMP(I))
           DENS(3*NRD+I)=YH1(ICOMP(I)) 
           DENS(4*NRD+I)=H*TP(1,ICOMP(I)) 
           DENS(5*NRD+I)=H2*YH2(ICOMP(I)) 
  69       CONTINUE
           XOLDD=XOLD
           HHH=H
C --- COMPUTE SOLUTION AND FIRST DERIVATIVE AT MID-POINT ----
           DO 473 J=2,KC
           DBLENJ=NJ(J)
           DO 473 L=J,2,-1
           FACTOR=(DBLENJ/NJ(L-1))**2-1.D0
           DO 473 I=1,NRD
           YSAFE(L-1,I)=YSAFE(L,I)+(YSAFE(L,I)-YSAFE(L-1,I))/FACTOR
           YSAFEP(L-1,I)=YSAFEP(L,I)+(YSAFEP(L,I)-YSAFEP(L-1,I))/FACTOR
 473       CONTINUE  
           DO 474 I=1,NRD
           DENS(6*NRD+I)=YSAFE(1,I)
 474       DENS(7*NRD+I)=H*YSAFEP(1,I)
C
           DO 180 KMI=2,KMIT 
C --- COMPUTE KMI-TH DERIVATIVE AT MID-POINT ----
             KBEG=KMI/2
             IF (KMI.EQ.2*KBEG) THEN
                IF (KMI.EQ.2) THEN
                    DO 570 I=1,NRD
 570                YSAFE(1,I)=(DZ(ICOMP(I))+FSAFE(1,I))/2.D0
                    KBEG=2
                END IF
                DO 375 KK=KBEG,KC
                FACNJ=0.5D0*(NJ(KK)/2.D0)**(KMI-2)  
                IPT=(KK-1)**2+KK+KMI/2-2
                DO 371 I=1,NRD
 371            YSAFE(KK,I)=(FSAFE(IPT,I)+FSAFE(IPT+1,I))*FACNJ
 375            CONTINUE 
             ELSE
                DO 575 KK=KBEG,KC
                FACNJ=(NJ(KK)/2.D0)**(KMI-2)  
                IPT=(KK-1)**2+KK+KBEG-1
                DO 571 I=1,NRD
 571            YSAFE(KK,I)=FSAFE(IPT,I)*FACNJ
 575            CONTINUE 
             END IF
C------          EXTRAPOLATION
             DO 373 J=KBEG+1,KC
             DBLENJ=NJ(J)
             DO 373 L=J,KBEG+1,-1
             FACTOR=(DBLENJ/NJ(L-1))**2-1.D0
             DO 373 I=1,NRD
             YSAFE(L-1,I)=YSAFE(L,I)+(YSAFE(L,I)-YSAFE(L-1,I))/FACTOR
 373         CONTINUE  
             KRN=(KMI+6)*NRD
             DO 374 I=1,NRD
 374         DENS(KRN+I)=YSAFE(KBEG,I)*H2
             IF (KMI.EQ.KMIT) GOTO 180
C ------          COMPUTE DIFFERENCES
             DO 66 KK=(KMI+1)/2,KC
             LBEG=(KK-1)**2+KMI-1 
             LEND=KK**2
             IF (KMI.EQ.2) LBEG=LBEG+1
             DO 64 L=LEND,LBEG,-1
             DO 64 I=1,NRD
  64         FSAFE(L,I)=FSAFE(L,I)-FSAFE(L-1,I)
             IF (KMI.EQ.2) THEN 
                L=LBEG-1
                DO 65 I=1,NRD
  65            FSAFE(L,I)=FSAFE(L,I)-DZ(ICOMP(I)) 
             END IF
  66         CONTINUE
 180      CONTINUE 
          CALL INTPO2(NRD,DENS,KMIT) 
C --- ESTIMATION OF INTERPOLATION ERROR  
          IF (IDERR.EQ.0.AND.KMIT.GE.1) THEN
            ERRINT=0.D0
            DO 187 I=1,NRD
 187        ERRINT=ERRINT+(DENS((KMIT+6)*NRD+I)/SCAL(ICOMP(I)))**2 
            ERRINT=SQRT(ERRINT/NRD)*ERRFAC(KMIT)
            HOPTDE=H/MAX((ERRINT)**(1.D0/(KMIT+6)),0.01D0) 
            IF (ERRINT.GT.10.D0) THEN  
             H=HOPTDE
             X=XOLD
             NREJCT=NREJCT+1
             REJECT=.TRUE.  
             GOTO 10 
            END IF 
          END IF
          DO 189 I=1,N
 189      DZ(I)=YH2(I)
       END IF
       DO 70 I=1,N
       Y(I)=T(1,I)
  70   YP(I)=TP(1,I)
       NACCPT=NACCPT+1  
       IF (IOUT.GE.1) THEN 
          CALL SOLOUT (NACCPT+1,XOLD,X,Y,YP,N,DENS,NCOM,ICOMP,NRD,
     &                 RPAR,IPAR,IRTRN)
          IF (IRTRN.LT.0) GOTO 120 
       END IF 
C --- COMPUTE OPTIMAL ORDER
       IF (KC.EQ.2) THEN
          KOPT=3  
          IF (REJECT) KOPT=2  
          GO TO 80
       END IF
       IF (KC.LE.K) THEN
          KOPT=KC 
          IF (W(KC-1).LT.W(KC)*FAC3) KOPT=KC-1  
          IF (W(KC).LT.W(KC-1)*FAC4) KOPT=MIN(KC+1,KM-1)
       ELSE 
          KOPT=KC-1
          IF (KC.GT.3.AND.W(KC-2).LT.W(KC-1)*FAC3) KOPT=KC-2
          IF (W(KC).LT.W(KOPT)*FAC4) KOPT=MIN(KC,KM-1)
       END IF
C --- AFTER A REJECTED STEP
  80   IF (REJECT) THEN 
          K=MIN(KOPT,KC)
          H=MIN(H,HH(K))
          REJECT=.FALSE.
          GO TO 10
       END IF
C --- COMPUTE STEPSIZE FOR NEXT STEP
       IF (KOPT.LE.KC) THEN
          H=HH(KOPT)
       ELSE 
          IF (KC.LT.K.AND.W(KC).LT.W(KC-1)*FAC4) THEN 
             H=HH(KC)*A(KOPT+1)/A(KC)
          ELSE
             H=HH(KC)*A(KOPT)/A(KC) 
          END IF  
       END IF
       K=KOPT
       IF (LAST) GOTO 110
       GO TO 10
C --- STEP IS REJECTED  
 100   CONTINUE
       K=MIN(K,KC)
       IF (K.GT.2.AND.W(K-1).LT.W(K)*FAC3) K=K-1
       NREJCT=NREJCT+1  
       H=HH(K)
       REJECT=.TRUE.
       LAST=.FALSE.
       GO TO 30
C --- SOLUTION EXIT
 110   CONTINUE 
       IDID=1
       RETURN
C --- FAIL EXIT
 120   WRITE (6,979) X,H
 979   FORMAT(' EXIT OF ODEX2 AT X=',D14.7,'   H=',D14.7)
       IDID=-1
       RETURN
       END  
C
       SUBROUTINE STOERM(J,X,Y,YP,H,HMAX,N,FCN,DY,YH1,YH2,DZ,T,TP,NJ,
     1   HH,W,ERR,FAC,A,SAFE1,UROUND,FAC1,FAC2,SAFE2,SCAL,SCALP,ATOV,
     2   SAFE3,REJECT,KM,RTOL,ATOL,ITOL,ERROLD,FSAFE,
     3   LFSAFE,IOUT,IPT,YSAFE,YSAFEP,ICOMP,NRD,RPAR,IPAR,NFCN)
C --- THIS SUBROUTINE COMPUTES THE J-TH LINE OF THE
C --- EXTRAPOLATION TABLE AND PROVIDES AN ESTIMATION  
C --- OF THE OPTIMAL STEPSIZE 
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       EXTERNAL FCN
       LOGICAL REJECT,ATOV
       DIMENSION Y(N),YP(N),DY(N),YH1(N),YH2(N),DZ(N),SCAL(N),SCALP(N) 
       DIMENSION T(KM,N),TP(KM,N),NJ(KM),HH(KM),W(KM),A(KM)
       DIMENSION RTOL(*),ATOL(*),RPAR(*),IPAR(*),ICOMP(NRD)
       DIMENSION FSAFE(LFSAFE,NRD),YSAFE(KM,NRD),YSAFEP(KM,NRD)
       HJ=H/NJ(J)
       HJ2=HJ*2.D0
C --- EULER STARTING STEP
       DO 30 I=1,N
       YH1(I)=Y(I)
  30   YH2(I)=YP(I)+HJ*DZ(I)
C --- EXPLICIT MIDPOINT (STOERMER) RULE  
       M=NJ(J)/2 
       DO 35 MM=1,M 
       IF (IOUT.GE.2.AND.MM.EQ.(M+1)/2) THEN
          DO 31 I=1,NRD
          YSAFEP(J,I)=YH2(ICOMP(I))
  31      YSAFE(J,I) =YH1(ICOMP(I))+HJ*YH2(ICOMP(I))
       END IF
       DO 33 I=1,N
  33   YH1(I)=YH1(I)+HJ2*YH2(I)
       CALL FCN(N,X+HJ2*MM,YH1,DY,RPAR,IPAR) 
       IF (IOUT.GE.2.AND.ABS(MM-(M+1)/2).LT.J) THEN 
          IPT=IPT+1
          DO 32 I=1,NRD
  32      FSAFE(IPT,I)=DY(ICOMP(I)) 
       END IF
       IF (MM.EQ.M) GOTO 37
       DO 34 I=1,N
  34   YH2(I)=YH2(I)+HJ2*DY(I)
  35   CONTINUE
C --- FINAL SMOOTHING STEP
  37   CONTINUE
       DO 40 I=1,N
       T(J,I)=YH1(I)
  40   TP(J,I)=YH2(I)+HJ*DY(I)
       NFCN=NFCN+M  
C --- POLYNOMIAL EXTRAPOLATION
       IF (J.EQ.1) RETURN
       DBLENJ=NJ(J)
       DO 60 L=J,2,-1
       FAC=(DBLENJ/NJ(L-1))**2-1.D0
       DO 60 I=1,N
       T(L-1,I)=T(L,I)+(T(L,I)-T(L-1,I))/FAC
       TP(L-1,I)=TP(L,I)+(TP(L,I)-TP(L-1,I))/FAC
  60   CONTINUE
       ERR=0.D0
C --- SCALING
       IF (ITOL.EQ.0) THEN
         DO 65 I=1,N
         SCAL(I)=ATOL(1)+RTOL(1)*MAX(ABS(Y(I)),ABS(T(1,I)))
  65     SCALP(I)=ATOL(1)+RTOL(1)*MAX(ABS(YP(I)),ABS(TP(1,I)))
       ELSE
         DO 66 I=1,N
         SCAL(I)=ATOL(I)+RTOL(I)*MAX(ABS(Y(I)),ABS(T(1,I)))
  66     SCALP(I)=ATOL(I+N)+RTOL(I+N)*MAX(ABS(YP(I)),ABS(TP(1,I)))
       END IF
       DO 67 I=1,N
  67   ERR=ERR+((T(1,I)-T(2,I))/SCAL(I))**2
     &        +((TP(1,I)-TP(2,I))/SCALP(I))**2  
       ERR=SQRT(ERR/(2*N))
       IF (ERR*UROUND.GE.1.D0) GOTO 79
       IF (J.GT.2.AND.ERR.GE.ERROLD) GOTO 79
       ERROLD=MAX(4*ERR,1.D0)
C --- COMPUTE OPTIMAL STEPSIZES
       EXPO=1.D0/(2*J-1)
       FACMIN=FAC1**EXPO
       FAC=MIN(FAC2/FACMIN,MAX(FACMIN,(ERR/SAFE1)**EXPO/SAFE2))
       FAC=1.D0/FAC
       HH(J)=MIN(H*FAC,HMAX)
       W(J)=A(J)/HH(J)  
       RETURN
  79   ATOV=.TRUE.
       H=H*SAFE3
       REJECT=.TRUE.
       RETURN
       END  
C
      SUBROUTINE INTPO2(N,Y,IMIT)
C --- COMPUTES THE COEFFICIENTS OF THE INTERPOLATION FORMULA
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(N*(IMIT+7)),A(0:40) 
C --- BEGIN WITH HERMITE INTERPOLATION (ORDER 5)
      DO 100 I=1,N
      Y0=Y(I)
      Y1=Y(3*N+I)
      YP0=Y(N+I)
      YP1=Y(4*N+I) 
      YPP0=Y(2*N+I)/2.D0
      YPP1=Y(5*N+I)/2.D0 
      YDIFF=Y1-Y0
      AH=YDIFF-YP0
      BH=YP1-YDIFF 
      CH=AH-YPP0
      DH=BH-AH
      EH=YPP1-BH
      FH=DH-CH
      GH=EH-DH
      Y(N+I)=YDIFF
      Y(2*N+I)=-AH
      Y(3*N+I)=-DH 
      Y(4*N+I)=GH 
      Y(5*N+I)=FH 
      IF (IMIT.LT.0) GOTO 100 
C --- COMPUTE THE DERIVATIVES OF HERMITE AT MIDPOINT
      ABH=AH+BH
      GFH=GH+FH
      GMF=GH-FH
      PH0=0.5D0*(Y0+Y1+0.25D0*(-ABH+0.25D0*GFH))
      PH1=YDIFF+0.25D0*(AH-BH+0.25D0*GMF)
      PH2=ABH-0.5D0*GFH
      PH3=6.D0*(BH-AH)-3.D0*GMF
      PH4=12.D0*GFH
      PH5=120.D0*GMF
C --- COMPUTE THE FURTHER COEFFICIENTS 
      IF (IMIT.LT.1) GOTO 20
      A(1)=64.D0*(Y(7*N+I)-PH1) 
      IF (IMIT.LT.3) GOTO 20 
      A(3)=64.D0*(Y(9*N+I)-PH3+A(1)*9.D0/8.D0)
      IF (IMIT.LT.5) GOTO 20 
      A(5)=64.D0*(Y(11*N+I)-PH5+A(3)*15.D0/4.D0-A(1)*90.D0)
      IF (IMIT.LT.7) GOTO 20 
      DO 10 IM=7,IMIT,2
      FC1=IM*(IM-1)*3.D0/16.D0
      FC2=FC1*(IM-2)*(IM-3)*4.D0
      FC3=IM*(IM-1)*(IM-2)*(IM-3)*(IM-4)*(IM-5)
  10  A(IM)=64.D0*(Y((IM+6)*N+I)+FC1*A(IM-2)-FC2*A(IM-4)+FC3*A(IM-6)) 
  20  CONTINUE 
      A(0)=64.D0*(Y(6*N+I)-PH0)
      IF (IMIT.LT.2) GOTO 60
      A(2)=64.D0*(Y(N*8+I)-PH2+A(0)*3.D0/8.D0)
      IF (IMIT.LT.4) GOTO 60
      A(4)=64.D0*(Y(N*10+I)-PH4+A(2)*9.D0/4.D0-A(0)*18.D0)
      IF (IMIT.LT.6) GOTO 60
      DO 30 IM=6,IMIT,2
      FC1=IM*(IM-1)*3.D0/16.D0
      FC2=FC1*(IM-2)*(IM-3)*4.D0
      FC3=IM*(IM-1)*(IM-2)*(IM-3)*(IM-4)*(IM-5)
  30  A(IM)=64.D0*(Y(N*(IM+6)+I)+A(IM-2)*FC1-A(IM-4)*FC2+A(IM-6)*FC3)
  60  CONTINUE 
      DO 70 IM=0,IMIT
  70  Y(N*(IM+6)+I)=A(IM)
 100  CONTINUE
      RETURN
      END 
C
      FUNCTION CONTX2(II,X,Y,NCON,ICOMP,N)
C ----------------------------------------------------------
C     THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT IN CONECTION
C     WITH THE OUTPUT-SUBROUTINE FOR ODEX2. IT PROVIDES AN
C     APPROXIMATION TO THE II-TH COMPONENT OF THE SOLUTION AT X.
C ----------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(NCON),ICOMP(N)
      COMMON /CONOD2/XOLD,H,IMIT
C ----- COMPUTE PLACE OF II-TH COMPONENT 
      I=0 
      DO 5 J=1,N 
      IF (ICOMP(J).EQ.II) I=J
   5  CONTINUE
      IF (I.EQ.0) THEN
         WRITE (6,*) ' NO DENSE OUTPUT AVAILABLE FOR COMP.',II 
         RETURN
      END IF  
C ----- COMPUTE THE INTERPOLATED VALUE 
      THETA=(X-XOLD)/H
      THETA1=1.D0-THETA
      PHTHET=Y(I)+THETA*(Y(N+I)+THETA1*(Y(2*N+I)+THETA*(Y(3*N+I)+
     &            THETA1*(Y(4*N+I)*THETA+Y(5*N+I)*THETA1))))
      IF (IMIT.LT.0) THEN
          CONTX2=PHTHET
          RETURN
      END IF
      THETAH=THETA-0.5D0
      CONTX2=Y(N*(IMIT+6)+I)
      DO 70 IM=IMIT,1,-1
  70  CONTX2=Y(N*(IM+5)+I)+CONTX2*THETAH/IM
      CONTX2=PHTHET+(THETA*THETA1)**3*CONTX2
      RETURN
      END