      SUBROUTINE RODAS(N,FCN,IFCN,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JAC ,IJAC,MLJAC,MUJAC,DFX,IDFX,
     &                  MAS ,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,LRCONT,IDID)
C ----------------------------------------------------------
C     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
C     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS  MY'=F(X,Y).
C     THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4  
C     (WITH STEP SIZE CONTROL).
C     C.F. SECTIONS IV.7  AND VI.3
C
C     AUTHORS: E. HAIRER AND G. WANNER
C              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
C              CH-1211 GENEVE 24, SWITZERLAND 
C              E-MAIL:  HAIRER@CGEUGE51.BITNET,  WANNER@CGEUGE51.BITNET
C     
C     THIS CODE IS PART OF THE BOOK:
C         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
C         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
C         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
C         SPRINGER-VERLAG (1990)               
C      
C     VERSION OF AUGUST 28, 1991
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)
C                    REAL*8 X,Y(N),F(N)
C                    F(1)=...   ETC.
C
C     IFCN        GIVES INFORMATION ON FCN:
C                    IFCN=0: F(X,Y) INDEPENDENT OF X (AUTONOMOUS)
C                    IFCN=1: F(X,Y) MAY DEPEND ON X (NON-AUTONOMOUS)
C
C     X           INITIAL X-VALUE
C
C     Y(N)        INITIAL VALUES FOR Y
C
C     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
C
C     H           INITIAL STEP SIZE GUESS;
C                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT, 
C                 H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD.
C                 THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY
C                 ADAPTS ITS STEP SIZE. STUDY THE CHOSEN VALUES FOR A FEW
C                 STEPS IN SUBROUTINE "SOLOUT", WHEN YOU ARE NOT SURE.
C                 (IF H=0.D0, THE CODE PUTS H=1.D-6).
C
C     RTOL,ATOL   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
C                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH 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                   ITOL=1: BOTH RTOL AND ATOL ARE VECTORS.
C                     THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW
C                     RTOL(I)*ABS(Y(I))+ATOL(I).
C
C     JAC         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
C                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y
C                 (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY
C                 A DUMMY SUBROUTINE IN THE CASE IJAC=0).
C                 FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM
C                    SUBROUTINE JAC(N,X,Y,DFY,LDFY)
C                    REAL*8 X,Y(N),DFY(LDFY,N)
C                    DFY(1,1)= ...
C                 LDFY, THE COLOMN-LENGTH OF THE ARRAY, IS
C                 FURNISHED BY THE CALLING PROGRAM.
C                 IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO
C                    BE FULL AND THE PARTIAL DERIVATIVES ARE
C                    STORED IN DFY AS
C                       DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J)
C                 ELSE, THE JACOBIAN IS TAKEN AS BANDED AND
C                    THE PARTIAL DERIVATIVES ARE STORED
C                    DIAGONAL-WISE AS
C                       DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).
C
C     IJAC        SWITCH FOR THE COMPUTATION OF THE JACOBIAN:
C                    IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE
C                       DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED.
C                    IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.
C
C     MLJAC       SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN:
C                    MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR
C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
C                    0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN 
C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
C                       THE MAIN DIAGONAL).
C
C     MUJAC       UPPER BANDWITH OF JACOBIAN  MATRIX (>= NUMBER OF NON-
C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
C                 NEED NOT BE DEFINED IF MLJAC=N.
C
C     DFX         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
C                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO X
C                 (THIS ROUTINE IS ONLY CALLED IF IDFX=1 AND IFCN=1;
C                 SUPPLY A DUMMY SUBROUTINE IN THE CASE IDFX=0 OR IFCN=0).
C                 OTHERWISE, THIS SUBROUTINE MUST HAVE THE FORM
C                    SUBROUTINE DFX(N,X,Y,FX)
C                    REAL*8 X,Y(N),FX(N)
C                    FX(1)= ...
C                
C     IDFX        SWITCH FOR THE COMPUTATION OF THE DF/DX:
C                    IDFX=0: DF/DX IS COMPUTED INTERNALLY BY FINITE
C                       DIFFERENCES, SUBROUTINE "DFX" IS NEVER CALLED.
C                    IDFX=1: DF/DX IS SUPPLIED BY SUBROUTINE DFX.
C
C     ----   MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS      -----
C     ----   FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): -
C
C     MAS         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS-
C                 MATRIX M.
C                 IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY
C                 MATRIX AND NEEDS NOT TO BE DEFINED;
C                 SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
C                 IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM
C                    SUBROUTINE MAS(N,AM,LMAS)
C                    REAL*8 AM(LMAS,N)
C                    AM(1,1)= ....
C                    IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED
C                    AS FULL MATRIX LIKE
C                         AM(I,J) = M(I,J)
C                    ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED
C                    DIAGONAL-WISE AS
C                         AM(I-J+MUMAS+1,J) = M(I,J).
C
C     IMAS       GIVES INFORMATION ON THE MASS-MATRIX:
C                    IMAS=0: M IS SUPPOSED TO BE THE IDENTITY
C                       MATRIX, MAS IS NEVER CALLED.
C                    IMAS=1: MASS-MATRIX  IS SUPPLIED.
C
C     MLMAS       SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX:
C                    MLMAS=N: THE FULL MATRIX CASE. THE LINEAR
C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
C                    0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE
C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
C                       THE MAIN DIAGONAL).
C                 MLMAS IS SUPPOSED TO BE .LE. MLJAC.
C
C     MUMAS       UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON-
C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
C                 NEED NOT BE DEFINED IF MLMAS=N.
C                 MUMAS IS SUPPOSED TO BE .LE. MUJAC.
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,N,IRTRN)
C                    REAL*8 X,Y(N)
C                    ....  
C                 SOLOUT FURNISHES THE SOLUTION "Y" 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, RODAS RETURNS TO THE CALLING PROGRAM.
C           
C          -----  CONTINUOUS OUTPUT: -----
C                 DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
C                 FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH
C                 THE FUNCTION
C                        >>>   CONTRO(I,S)   <<<
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        GIVES INFORMATION ON THE SUBROUTINE SOLOUT:
C                    IOUT=0: SUBROUTINE IS NEVER CALLED
C                    IOUT=1: SUBROUTINE IS USED FOR OUTPUT
C
C     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
C                 SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES.
C                 "LWORK" MUST BE AT LEAST
C                             N*(LJAC+LMAS+LE1+10)+4
C                 WHERE
C                    LJAC=N              IF MLJAC=N (FULL JACOBIAN)
C                    LJAC=MLJAC+MUJAC+1  IF MLJAC<N (BANDED JAC.)
C                 AND                  
C                    LMAS=0              IF IMAS=0
C                    LMAS=N              IF IMAS=1 AND MLMAS=N (FULL)
C                    LMAS=MLMAS+MUMAS+1  IF MLMAS<N (BANDED MASS-M.)
C                 AND
C                    LE1=N               IF MLJAC=N (FULL JACOBIAN)
C                    LE1=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.).
C                 IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE
C                 MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM
C                 STORAGE REQUIREMENT IS 
C                             LWORK = 2*N*N+10*N+4.
C
C     LWORK       DECLARED LENGHT OF ARRAY "WORK".
C
C     IWORK       INTEGER WORKING SPACE OF LENGHT "LIWORK".
C                 "LIWORK" MUST BE AT LEAST N+2.
C
C     LIWORK      DECLARED LENGHT OF ARRAY "IWORK".
C
C     LRCONT      DECLARED LENGTH OF COMMON BLOCK
C                  >>>  COMMON /CONROS/ICONT(3),RCONT(LRCONT)  <<<
C                 WHICH MUST BE DECLARED IN THE CALLING PROGRAM.
C                 "LRCONT" MUST BE AT LEAST
C                             4*N+2 .
C                 THIS IS USED FOR STORING THE COEFFICIENTS OF THE
C                 CONTINUOUS SOLUTION AND MAKES THE CALLING LIST FOR THE
C                 FUNCTION "CONTD5" AS SIMPLE AS POSSIBLE.
C
C ----------------------------------------------------------------------
C 
C     SOPHISTICATED SETTING OF PARAMETERS
C     -----------------------------------
C              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK 
C              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(4)
C              AS WELL AS IWORK(1),IWORK(2) DIFFERENT FROM ZERO.
C              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
C
C    IWORK(1)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
C              THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000.
C
C    IWORK(2)  SWITCH FOR THE CHOICE OF THE COEFFICIENTS
C              IF IWORK(2).EQ.1  METHOD WITH ?????
C              IF IWORK(2).EQ.2  METHOD WITH ?????
C              THE DEFAULT VALUE (FOR IWORK(2)=0) IS IWORK(2)=1.
C
C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
C
C    WORK(2)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
C
C    WORK(3), WORK(4)   PARAMETERS FOR STEP SIZE SELECTION
C              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION
C                 WORK(3) <= HNEW(J)/HOLD <= WORK(4)
C              DEFAULT VALUES: WORK(3)=0.2D0, WORK(4)=6.D0
C
C-----------------------------------------------------------------------
C
C     OUTPUT PARAMETERS 
C     ----------------- 
C     X           X-VALUE WHERE THE SOLUTION IS COMPUTED
C                 (AFTER SUCCESSFUL RETURN X=XEND)
C
C     Y(N)        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 --------------------------------------------------------- 
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C          DECLARATIONS 
C *** *** *** *** *** *** *** *** *** *** *** *** ***
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION Y(N),ATOL(1),RTOL(1),WORK(LWORK),IWORK(LIWORK)
      LOGICAL AUTNMS,IMPLCT,JBAND,ARRET
      EXTERNAL FCN,JAC,DFX,MAS,SOLOUT
      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
C -----------------------------------------------------
C --- COMMON STAT CAN BE USED FOR STATISTICS
C ---    NFCN      NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL
C                  EVALUATION OF THE JACOBIAN ARE NOT COUNTED)  
C ---    NJAC      NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
C                  OR NUMERICALLY)
C ---    NSTEP     NUMBER OF COMPUTED STEPS
C ---    NACCPT    NUMBER OF ACCEPTED STEPS
C ---    NREJCT    NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP
C                  HAS BEEN ACCEPTED)
C ---    NDEC      NUMBER OF LU-DECOMPOSITIONS (N-DIMENSIONAL MATRIX)
C ---    NSOL      NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS
C -----------------------------------------------------------
C *** *** *** *** *** *** ***
C        SETTING THE PARAMETERS 
C *** *** *** *** *** *** ***
      NFCN=0
      NACCPT=0
      NREJCT=0
      NSTEP=0
      NJAC=0
      NDEC=0
      NSOL=0
      ARRET=.FALSE.
C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
      IF(IWORK(1).EQ.0)THEN
         NMAX=100000
      ELSE
         NMAX=IWORK(1)
         IF(NMAX.LE.0)THEN
            WRITE(6,*)' WRONG INPUT IWORK(1)=',IWORK(1)
            ARRET=.TRUE.
         END IF
      END IF
C -------- METH   COEFFICIENTS OF THE METHOD
      IF(IWORK(2).EQ.0)THEN
         METH=1
      ELSE
         METH=IWORK(2)
         IF(METH.LE.0.OR.METH.GE.3)THEN
            WRITE(6,*)' CURIOUS INPUT IWORK(2)=',IWORK(2)
            ARRET=.TRUE.
         END IF
      END IF
C -------- UROUND   SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0  
      IF(WORK(1).EQ.0.D0)THEN
         UROUND=1.D-16
      ELSE
         UROUND=WORK(1)
         IF(UROUND.LE.1.D-14.OR.UROUND.GE.1.D0)THEN
            WRITE(6,*)' COEFFICIENTS HAVE 16 DIGITS, UROUND=',WORK(1)
            ARRET=.TRUE.
         END IF
      END IF
C -------- MAXIMAL STEP SIZE
      IF(WORK(2).EQ.0.D0)THEN
         HMAX=XEND-X
      ELSE
         HMAX=WORK(2)
      END IF
C -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
      IF(WORK(3).EQ.0.D0)THEN
         FAC1=5.D0
      ELSE
         FAC1=1.D0/WORK(3)
      END IF
      IF(WORK(4).EQ.0.D0)THEN
         FAC2=1.D0/6.0D0
      ELSE
         FAC2=1.D0/WORK(4)
      END IF
C --------- CHECK IF TOLERANCES ARE O.K.
      IF (ITOL.EQ.0) THEN
          IF (ATOL(1).LE.0.D0.OR.RTOL(1).LE.10.D0*UROUND) THEN
              WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
              ARRET=.TRUE.
          END IF
      ELSE
          DO 15 I=1,N
          IF (ATOL(I).LE.0.D0.OR.RTOL(I).LE.10.D0*UROUND) THEN
              WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
              ARRET=.TRUE.
          END IF
  15      CONTINUE
      END IF
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C         COMPUTATION OF ARRAY ENTRIES
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ?
      AUTNMS=IFCN.EQ.0
      IMPLCT=IMAS.NE.0
      JBAND=MLJAC.NE.N
      ARRET=.FALSE.
C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
C -- JACOBIAN AND MATRIX E
      IF(JBAND)THEN
         LDJAC=MLJAC+MUJAC+1
         LDE=MLJAC+LDJAC
      ELSE
         LDJAC=N
         LDE=N
      END IF
C -- MASS MATRIX
      IF (IMPLCT) THEN
          IF (MLMAS.NE.N) THEN
              LDMAS=MLMAS+MUMAS+1
          ELSE
              LDMAS=N
          END IF
C ------ BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF "JAC"
          IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN
              WRITE (6,*) 'BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF
     & "JAC"'
            ARRET=.TRUE.
          END IF
      ELSE
          LDMAS=0
      END IF
      LDMAS2=MAX(1,LDMAS)
C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
      IEYNEW=5
      IEDY1=IEYNEW+N
      IEDY=IEDY1+N
      IEAK1=IEDY+N
      IEAK2=IEAK1+N
      IEAK3=IEAK2+N
      IEAK4=IEAK3+N
      IEAK5=IEAK4+N
      IEAK6=IEAK5+N
      IEFX =IEAK6+N
      IEJAC=IEFX +N
      IEMAS=IEJAC+N*LDJAC
      IEE  =IEMAS+N*LDMAS
C ------ TOTAL STORAGE REQUIREMENT -----------
      ISTORE=IEE+N*LDE-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 -----
      IEIP=3
      ISTORE=IEIP+N-1
      IF(ISTORE.GT.LIWORK)THEN
         WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
         ARRET=.TRUE.
      END IF
C --------- REQUIREMENT FOR RCONT ---------------
      IF(4*N+2.GT.LRCONT)THEN
         WRITE(6,*)' INSUFF. STORAGE FOR RCONT, MIN. LRCONT=',4*N+2
         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 ------------
      CALL ROSCOR(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,IJAC,
     &   MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
     &   NMAX,UROUND,METH,FAC1,FAC2,AUTNMS,IMPLCT,JBAND,LDJAC,
     &   LDE,LDMAS2,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY),WORK(IEAK1),
     &   WORK(IEAK2),WORK(IEAK3),WORK(IEAK4),WORK(IEAK5),WORK(IEAK6),
     &   WORK(IEFX),WORK(IEJAC),WORK(IEE),WORK(IEMAS),IWORK(IEIP))
C ----------- RETURN -----------
      RETURN
      END
C
C
C
C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
C
      SUBROUTINE ROSCOR(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,
     &  IJAC,MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
     &  NMAX,UROUND,METH,FAC1,FAC2,AUTNMS,IMPLCT,BANDED,LDJAC,LDE,
     &  LDMAS,YNEW,DY1,DY,AK1,AK2,AK3,AK4,AK5,AK6,FX,FJAC,E,FMAS,IP)
C ----------------------------------------------------------
C     CORE INTEGRATOR FOR RODAS
C     PARAMETERS SAME AS IN RODAS WITH WORKSPACE ADDED 
C ---------------------------------------------------------- 
C         DECLARATIONS 
C ---------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Y(N),YNEW(N),DY1(N),DY(N),AK1(N),
     *  AK2(N),AK3(N),AK4(N),AK5(N),AK6(N),FX(N),
     *  FJAC(LDJAC,N),E(LDE,N),FMAS(LDMAS,N),ATOL(1),RTOL(1)
      INTEGER IP(N)
      LOGICAL REJECT,AUTNMS,IMPLCT,BANDED
      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
      COMMON /CONROS/NN,NN2,NN3,XOLD,HOUT,CONT(1)
C *** *** *** *** *** *** ***
C  INITIALISATIONS
C *** *** *** *** *** *** ***    
      NN=N 
      NN2=2*N
      NN3=3*N
C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
      IF(IMPLCT)CALL MAS(N,FMAS,LDMAS)
C ------- FORWAED OR BACKWARD INTEGRATION ---- 
      POSNEG=DSIGN(1.D0,XEND-X)
C ------ SET THE PARAMETERS OF THE METHOD -----
      IF(METH.EQ.1)CALL ROCOE1(A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,
     &  C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,
     &  C62,C63,C64,C65,GAMMA,C2,C3,C4,D1,D2,D3,D4,
     &  D21,D22,D23,D24,D25,D31,D32,D33,D34,D35)
      IF(METH.EQ.2)CALL ROCOE2(A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,
     &  C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,
     &  C62,C63,C64,C65,GAMMA,C2,C3,C4,D1,D2,D3,D4,
     &  D21,D22,D23,D24,D25,D31,D32,D33,D34,D35)
C --- INITIAL PREPARATIONS
      HMAXN=DMIN1(DABS(HMAX),DABS(XEND-X))
      H=DMIN1(DMAX1(1.D-10,DABS(H)),HMAXN) 
      H=DSIGN(H,POSNEG) 
      REJECT=.FALSE.
      NSING=0
      IRTRN=1
C -------- PREPARE BAND-WIDTHS --------
      MBDIAG=MUMAS+1
      IF (BANDED) THEN
          MLE=MLJAC
          MUE=MUJAC
          MBJAC=MLJAC+MUJAC+1
          MBB=MLMAS+MUMAS+1
          MDIAG=MLE+MUE+1
          MDIFF=MLE+MUE-MUMAS
      END IF
      IF (IOUT.NE.0) THEN 
          XOLD=X
          IRTRN=1
          HOUT=H
          CALL SOLOUT(NACCPT+1,XOLD,X,Y,N,IRTRN)
          IF (IRTRN.LT.0) GOTO 79
      END IF
C --- BASIC INTEGRATION STEP  
   1  IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79
      IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN
          H=HOPT
          IDID=1
          RETURN
      END IF
      HOPT=H
      IF ((X+H-XEND)*POSNEG.GT.0.D0) H=XEND-X
C *** *** *** *** *** *** ***
C  COMPUTATION OF THE JACOBIAN
C *** *** *** *** *** *** ***
      CALL FCN(N,X,Y,DY1)
      NFCN=NFCN+1
      NJAC=NJAC+1
      IF (IJAC.EQ.0) THEN
C --- COMPUTE JACOBIAN MATRIX NUMERICALLY
          IF (BANDED) THEN
C --- JACOBIAN IS BANDED
              MUJACP=MUJAC+1
              MD=MIN(MBJAC,N)
              DO 16 K=1,MD
              J=K
 12           AK2(J)=Y(J)
              AK3(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J))))
              Y(J)=Y(J)+AK3(J)
              J=J+MD
              IF (J.LE.N) GOTO 12 
              CALL FCN(N,X,Y,AK1)
              J=K
              LBEG=MAX(1,J-MUJAC)
 14           LEND=MIN(N,J+MLJAC)
              Y(J)=AK2(J)
              MUJACJ=MUJACP-J
              DO 15 L=LBEG,LEND
 15           FJAC(L+MUJACJ,J)=(AK1(L)-DY1(L))/AK3(J) 
              J=J+MD
              LBEG=LEND+1
              IF (J.LE.N) GOTO 14
 16           CONTINUE
          ELSE
C --- JACOBIAN IS FULL
              DO 18 I=1,N
              YSAFE=Y(I)
              DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE)))
              Y(I)=YSAFE+DELT
              CALL FCN(N,X,Y,AK1)
              DO 17 J=1,N
  17          FJAC(J,I)=(AK1(J)-DY1(J))/DELT
  18          Y(I)=YSAFE
              MLJAC=N
          END IF
      ELSE
C --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
          CALL JAC(N,X,Y,FJAC,LDJAC)
      END IF
      IF (.NOT.AUTNMS) THEN
          IF (IDFX.EQ.0) THEN
C --- COMPUTE NUMERICALLY THE DERIVATIVE WITH RESPECT TO X
              DELT=DSQRT(UROUND*MAX(1.D-5,ABS(X)))
              XDELT=X+DELT
              CALL FCN(N,XDELT,Y,AK1)
              DO 19 J=1,N
  19          FX(J)=(AK1(J)-DY1(J))/DELT
          ELSE
C --- COMPUTE ANALYTICALLY THE DERIVATIVE WITH RESPECT TO X
              CALL DFX(N,X,Y,FX)
          END IF
      END IF
   2  CONTINUE
C *** *** *** *** *** *** ***
C  COMPUTE THE STAGES
C *** *** *** *** *** *** ***
      NDEC=NDEC+1
C --- PREPARE FOR THE COMPUTATION OF THE 6 STAGES
      HC21=C21/H
      HC31=C31/H
      HC32=C32/H
      HC41=C41/H
      HC42=C42/H
      HC43=C43/H
      HC51=C51/H
      HC52=C52/H
      HC53=C53/H
      HC54=C54/H
      HC61=C61/H
      HC62=C62/H
      HC63=C63/H
      HC64=C64/H
      HC65=C65/H
      FAC=1.D0/(H*GAMMA)
      IF (IMPLCT) THEN
          IF (BANDED) THEN
C + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
C --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
C ---   1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM
C ---   2) THE MATRIX B AND THE JACOBIAN OF F ARE BANDED
C + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
C --- THE MATRIX E (B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX)
              DO 101 J=1,N
              I1=MAX0(1,MUJAC+2-J)
              I2=MIN0(MBJAC,N+MUJAC+1-J)
              DO 101 I=I1,I2
  101         E(I+MLE,J)=-FJAC(I,J)
              DO 102 J=1,N
              I1=MAX0(1,MUMAS+2-J)
              I2=MIN0(MBB,N+MUMAS+1-J)
              DO 102 I=I1,I2
              IB=I+MDIFF
  102         E(IB,J)=E(IB,J)+FAC*FMAS(I,J)
              CALL DECB(N,LDE,E,MLE,MUE,IP,INFO)
              IF (INFO.NE.0) GOTO 80
              IF (AUTNMS) THEN
C --- THE DIFFERENTIAL EQUATION IS AUTONOMOUS
                  DO 103 I=1,N
  103             AK1(I)=DY1(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK1,IP)
                  DO 110 I=1,N
  110             YNEW(I)=Y(I)+A21*AK1(I) 
                      CALL FCN(N,X,YNEW,DY)
                  DO 111 I=1,N
  111             YNEW(I)=HC21*AK1(I)
                  DO 114 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 113 J=J1,J2
  113             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  114             AK2(I)=SUM+DY(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK2,IP)
                  DO 120 I=1,N
  120             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                      CALL FCN(N,X,YNEW,DY)
                  DO 121 I=1,N
  121             YNEW(I)=HC31*AK1(I)+HC32*AK2(I)
                  DO 124 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 123 J=J1,J2
  123             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  124             AK3(I)=SUM+DY(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK3,IP)
                  DO 130 I=1,N
  130             YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)  
                      CALL FCN(N,X,YNEW,DY)
                  DO 131 I=1,N
  131             YNEW(I)=HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 
                  DO 134 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 133 J=J1,J2
  133             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  134             AK4(I)=SUM+DY(I)
             ELSE
C --- THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS
                  HD1=H*D1
                  HD2=H*D2
                  HD3=H*D3
                  HD4=H*D4
                  DO 203 I=1,N
  203             AK1(I)=DY1(I)+HD1*FX(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK1,IP)
                 DO 210 I=1,N
  210             YNEW(I)=Y(I)+A21*AK1(I) 
                      CALL FCN(N,X+C2*H,YNEW,DY)
                  DO 211 I=1,N
  211             YNEW(I)=HC21*AK1(I)
                  DO 214 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 213 J=J1,J2
  213             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  214             AK2(I)=SUM+DY(I)+HD2*FX(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK2,IP)
                  DO 220 I=1,N
  220             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                      CALL FCN(N,X+C3*H,YNEW,DY)
                  DO 221 I=1,N
  221             YNEW(I)=HC31*AK1(I)+HC32*AK2(I)
                  DO 224 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 223 J=J1,J2
  223             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  224             AK3(I)=SUM+DY(I)+HD3*FX(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK3,IP)
                  DO 230 I=1,N
  230             YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)  
                      CALL FCN(N,X+C4*H,YNEW,DY)
                  DO 231 I=1,N
  231             YNEW(I)=HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 
                  DO 234 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 233 J=J1,J2
  233             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  234             AK4(I)=SUM+DY(I)+HD4*FX(I)
              END IF
C --- THE DIFFERENTIAL EQUATION IS AUTONOMOUS OR NON-AUTONOMOUS
                  CALL SOLB(N,LDE,E,MLE,MUE,AK4,IP)
              DO 240 I=1,N
  240         YNEW(I)=Y(I)+A51*AK1(I)+A52*AK2(I)+A53*AK3(I)
     &                +A54*AK4(I)  
              CALL FCN(N,X+H,YNEW,DY)
              DO 241 I=1,N
  241         AK6(I)=HC51*AK1(I)+HC52*AK2(I)+HC53*AK3(I) 
     &              +HC54*AK4(I) 
              DO 244 I=1,N
              SUM=0.D0
              J1=MAX0(1,I-MLMAS)
              J2=MIN0(N,I+MUMAS)
              DO 243 J=J1,J2
  243         SUM=SUM+FMAS(I-J+MBDIAG,J)*AK6(J)
  244         AK5(I)=SUM+DY(I)
                  CALL SOLB(N,LDE,E,MLE,MUE,AK5,IP)
C ------------ EMBEDDED SOLUTION ---------------
              DO 250 I=1,N
  250         YNEW(I)=YNEW(I)+AK5(I)  
              CALL FCN(N,X+H,YNEW,DY)
              DO 251 I=1,N
  251         AK5(I)=HC61*AK1(I)+HC62*AK2(I)+HC63*AK3(I)
     &              +HC64*AK4(I)+HC65*AK5(I) 
              DO 254 I=1,N
              SUM=0.D0
              J1=MAX0(1,I-MLMAS)
              J2=MIN0(N,I+MUMAS)
              DO 253 J=J1,J2
  253         SUM=SUM+FMAS(I-J+MBDIAG,J)*AK5(J)
  254         AK6(I)=SUM+DY(I)
                  CALL SOLB(N,LDE,E,MLE,MUE,AK6,IP)
C ------------ NEW SOLUTION ---------------
              DO 255 I=1,N
  255         YNEW(I)=YNEW(I)+AK6(I)  
          ELSE
              IF (MLMAS.NE.N) THEN
C + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
C --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
C ---   1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM
C ---   2) THE MATRIX B IS BANDED BUT THE JACOBIAN OF F IS NOT
C + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
C --- THE MATRIX E (B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX)
                  DO 301 J=1,N
                  DO 301 I=1,N
  301             E(I,J)=-FJAC(I,J)
                  DO 302 J=1,N
                  I1=MAX0(1,J-MUMAS)
                  I2=MIN0(N,J+MLMAS)
                  DO 302 I=I1,I2
  302             E(I,J)=E(I,J)+FAC*FMAS(I-J+MBDIAG,J)
                  CALL DEC(N,LDE,E,IP,INFO)
                  IF (INFO.NE.0) GOTO 80
                  IF (AUTNMS) THEN
C --- THE DIFFERENTIAL EQUATION IS AUTONOMOUS
                  DO 303 I=1,N
  303             AK1(I)=DY1(I)
                      CALL SOL(N,LDE,E,AK1,IP)
                  DO 310 I=1,N
  310             YNEW(I)=Y(I)+A21*AK1(I) 
                      CALL FCN(N,X,YNEW,DY)
                  DO 311 I=1,N
  311             YNEW(I)=HC21*AK1(I)
                  DO 314 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 313 J=J1,J2
  313             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  314             AK2(I)=SUM+DY(I)
                      CALL SOL(N,LDE,E,AK2,IP)
                  DO 320 I=1,N
  320             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                      CALL FCN(N,X,YNEW,DY)
                  DO 321 I=1,N
  321             YNEW(I)=HC31*AK1(I)+HC32*AK2(I)
                  DO 324 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 323 J=J1,J2
  323             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  324             AK3(I)=SUM+DY(I)
                      CALL SOL(N,LDE,E,AK3,IP)
                  DO 330 I=1,N
  330             YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)  
                      CALL FCN(N,X,YNEW,DY)
                  DO 331 I=1,N
  331             YNEW(I)=HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 
                  DO 334 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 333 J=J1,J2
  333             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  334             AK4(I)=SUM+DY(I)
                ELSE
C --- THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS
                  HD1=H*D1
                  HD2=H*D2
                  HD3=H*D3
                  HD4=H*D4
                  DO 304 I=1,N
  304             AK1(I)=DY1(I)+HD1*FX(I)
                      CALL SOL(N,LDE,E,AK1,IP)
                  DO 305 I=1,N
  305             YNEW(I)=Y(I)+A21*AK1(I) 
                      CALL FCN(N,X+C2*H,YNEW,DY)
                  DO 306 I=1,N
  306             YNEW(I)=HC21*AK1(I)
                  DO 308 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 307 J=J1,J2
  307             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  308             AK2(I)=SUM+DY(I)+HD2*FX(I)
                      CALL SOL(N,LDE,E,AK2,IP)
                  DO 309 I=1,N
  309             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                      CALL FCN(N,X+C3*H,YNEW,DY)
                  DO 355 I=1,N
  355             YNEW(I)=HC31*AK1(I)+HC32*AK2(I)
                  DO 357 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 356 J=J1,J2
  356             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  357             AK3(I)=SUM+DY(I)+HD3*FX(I)
                      CALL SOL(N,LDE,E,AK3,IP)
                  DO 360 I=1,N
  360             YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)  
                      CALL FCN(N,X+C4*H,YNEW,DY)
                  DO 361 I=1,N
  361             YNEW(I)=HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 
                  DO 364 I=1,N
                  SUM=0.D0
                  J1=MAX0(1,I-MLMAS)
                  J2=MIN0(N,I+MUMAS)
                  DO 363 J=J1,J2
  363             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
  364             AK4(I)=SUM+DY(I)+HD4*FX(I)
                END IF
C --- THE DIFFERENTIAL EQUATION IS AUTONOMOUS OR NON-AUTONOMOUS
                    CALL SOL(N,LDE,E,AK4,IP)
                DO 370 I=1,N
  370           YNEW(I)=Y(I)+A51*AK1(I)+A52*AK2(I)+A53*AK3(I)
     &                +A54*AK4(I)  
                CALL FCN(N,X+H,YNEW,DY)
                DO 371 I=1,N
  371           AK6(I)=HC51*AK1(I)+HC52*AK2(I)+HC53*AK3(I) 
     &                +HC54*AK4(I) 
                DO 374 I=1,N
                SUM=0.D0
                J1=MAX0(1,I-MLMAS)
                J2=MIN0(N,I+MUMAS)
                DO 373 J=J1,J2
  373           SUM=SUM+FMAS(I-J+MBDIAG,J)*AK6(J)
  374           AK5(I)=SUM+DY(I)
                    CALL SOL(N,LDE,E,AK5,IP)
C ------------ EMBEDDED SOLUTION ---------------
                DO 380 I=1,N
  380           YNEW(I)=YNEW(I)+AK5(I)  
                CALL FCN(N,X+H,YNEW,DY)
                DO 381 I=1,N
  381           AK5(I)=HC61*AK1(I)+HC62*AK2(I)+HC63*AK3(I)
     &                +HC64*AK4(I)+HC65*AK5(I) 
                DO 384 I=1,N
                SUM=0.D0
                J1=MAX0(1,I-MLMAS)
                J2=MIN0(N,I+MUMAS)
                DO 383 J=J1,J2
  383           SUM=SUM+FMAS(I-J+MBDIAG,J)*AK5(J)
  384           AK6(I)=SUM+DY(I)
                    CALL SOL(N,LDE,E,AK6,IP)
C ------------ NEW SOLUTION ---------------
              DO 385 I=1,N
  385         YNEW(I)=YNEW(I)+AK6(I)  
               ELSE
C + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
C --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
C ---   1) THE DIFFERENTIAL EQUATION IS IN IMPLICIT FORM
C ---   2) THE MATRIX B IS NOT BANDED 
C + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
C --- THE MATRIX E (B IS A FULL MATRIX, JACOBIAN A FULL OR BANDED MATRIX)
                  IF (MLJAC.EQ.N) THEN
                      DO 401 J=1,N
                      DO 401 I=1,N
  401                 E(I,J)=FMAS(I,J)*FAC-FJAC(I,J)
                  ELSE
                      MADD=MUJAC+1
                      DO 405 J=1,N
                      DO 405 I=1,N
  405                 E(I,J)=FMAS(I,J)*FAC
                      DO 406 J=1,N
                      I1=MAX0(1,J-MUJAC)
                      I2=MIN0(N,J+MLJAC)
                      DO 406 I=I1,I2
  406                 E(I,J)=E(I,J)-FJAC(I-J+MADD,J)
                  END IF
                  CALL DEC(N,LDE,E,IP,INFO)
                  IF (INFO.NE.0) GOTO 80
                  IF (AUTNMS) THEN
C --- THE DIFFERENTIAL EQUATION IS AUTONOMOUS
                      DO 403 I=1,N
  403                 AK1(I)=DY1(I)
                          CALL SOL(N,LDE,E,AK1,IP)
                      DO 410 I=1,N
  410                 YNEW(I)=Y(I)+A21*AK1(I) 
                          CALL FCN(N,X,YNEW,DY)
                      DO 411 I=1,N
  411                 YNEW(I)=HC21*AK1(I)
                      DO 414 I=1,N
                      SUM=0.D0
                      DO 413 J=1,N
  413                 SUM=SUM+FMAS(I,J)*YNEW(J)
  414                 AK2(I)=SUM+DY(I)
                      CALL SOL(N,LDE,E,AK2,IP)
                      DO 420 I=1,N
  420                 YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                          CALL FCN(N,X,YNEW,DY)
                      DO 421 I=1,N
  421                 YNEW(I)=HC31*AK1(I)+HC32*AK2(I)
                      DO 424 I=1,N
                      SUM=0.D0
                      DO 423 J=1,N
  423                 SUM=SUM+FMAS(I,J)*YNEW(J)
  424                 AK3(I)=SUM+DY(I)
                          CALL SOL(N,LDE,E,AK3,IP)
                      DO 430 I=1,N
  430                 YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)  
                          CALL FCN(N,X,YNEW,DY)
                      DO 431 I=1,N
  431                 YNEW(I)=HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 
                      DO 434 I=1,N
                      SUM=0.D0
                      DO 433 J=1,N
  433                 SUM=SUM+FMAS(I,J)*YNEW(J)
  434                 AK4(I)=SUM+DY(I)
                  ELSE
C --- THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS
                      HD1=H*D1
                      HD2=H*D2
                      HD3=H*D3
                      HD4=H*D4
                      DO 503 I=1,N
  503                 AK1(I)=DY1(I)+HD1*FX(I)
                          CALL SOL(N,LDE,E,AK1,IP)
                      DO 510 I=1,N
  510                 YNEW(I)=Y(I)+A21*AK1(I) 
                          CALL FCN(N,X+C2*H,YNEW,DY)
                      DO 511 I=1,N
  511                 YNEW(I)=HC21*AK1(I)
                      DO 514 I=1,N
                      SUM=0.D0
                      DO 513 J=1,N
  513                 SUM=SUM+FMAS(I,J)*YNEW(J)
  514                 AK2(I)=SUM+DY(I)+HD2*FX(I)
                      CALL SOL(N,LDE,E,AK2,IP)
                      DO 520 I=1,N
  520                 YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                          CALL FCN(N,X+C3*H,YNEW,DY)
                      DO 521 I=1,N
  521                 YNEW(I)=HC31*AK1(I)+HC32*AK2(I)
                      DO 524 I=1,N
                      SUM=0.D0
                      DO 523 J=1,N
  523                 SUM=SUM+FMAS(I,J)*YNEW(J)
  524                 AK3(I)=SUM+DY(I)+HD3*FX(I)
                          CALL SOL(N,LDE,E,AK3,IP)
                      DO 530 I=1,N
  530                 YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)  
                          CALL FCN(N,X+C4*H,YNEW,DY)
                      DO 531 I=1,N
  531                 YNEW(I)=HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 
                      DO 534 I=1,N
                      SUM=0.D0
                      DO 533 J=1,N
  533                 SUM=SUM+FMAS(I,J)*YNEW(J)
  534                 AK4(I)=SUM+DY(I)+HD4*FX(I)
                  END IF
C --- THE DIFFERENTIAL EQUATION IS AUTONOMOUS OR NON-AUTONOMOUS
                      CALL SOL(N,LDE,E,AK4,IP)
                  DO 440 I=1,N
  440             YNEW(I)=Y(I)+A51*AK1(I)+A52*AK2(I)+A53*AK3(I)
     &                  +A54*AK4(I)  
                  CALL FCN(N,X+H,YNEW,DY)
                  DO 441 I=1,N
  441             AK6(I)=HC51*AK1(I)+HC52*AK2(I)+HC53*AK3(I) 
     &                  +HC54*AK4(I) 
                  DO 444 I=1,N
                  SUM=0.D0
                  DO 443 J=1,N
  443             SUM=SUM+FMAS(I,J)*AK6(J)
  444             AK5(I)=SUM+DY(I)
                      CALL SOL(N,LDE,E,AK5,IP)
C ------------ EMBEDDED SOLUTION ---------------
                  DO 450 I=1,N
  450             YNEW(I)=YNEW(I)+AK5(I)  
                  CALL FCN(N,X+H,YNEW,DY)
                  DO 451 I=1,N
  451             AK5(I)=HC61*AK1(I)+HC62*AK2(I)+HC63*AK3(I)
     &                  +HC64*AK4(I)+HC65*AK5(I) 
                  DO 454 I=1,N
                  SUM=0.D0
                  DO 453 J=1,N
  453             SUM=SUM+FMAS(I,J)*AK5(J)
  454             AK6(I)=SUM+DY(I)
                      CALL SOL(N,LDE,E,AK6,IP)
C ------------ NEW SOLUTION ---------------
                  DO 455 I=1,N
  455             YNEW(I)=YNEW(I)+AK6(I)  
              END IF
          END IF
      ELSE
          IF (BANDED) THEN
C + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
C --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
C ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
C ---   2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX
C + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
C --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX)
              MLE=MLJAC
              MUE=MUJAC
              MBJAC=MLJAC+MUJAC+1
              MDIAG=MLE+MUE+1
              DO 601 J=1,N
              I1=MAX0(1,MUJAC+2-J)
              I2=MIN0(MBJAC,N+MUJAC+1-J)
              DO 600 I=I1,I2
  600         E(I+MLE,J)=-FJAC(I,J)
  601         E(MDIAG,J)=E(MDIAG,J)+FAC
              CALL DECB(N,LDE,E,MLE,MUE,IP,INFO)
              IF (INFO.NE.0) GOTO 80
              IF (AUTNMS) THEN
C --- THE DIFFERENTIAL EQUATION IS AUTONOMOUS
                  DO 603 I=1,N
  603             AK1(I)=DY1(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK1,IP)
                  DO 610 I=1,N
  610             YNEW(I)=Y(I)+A21*AK1(I) 
                  CALL FCN(N,X,YNEW,DY)
                  DO 611 I=1,N
  611             AK2(I)=DY(I)+HC21*AK1(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK2,IP)
                  DO 620 I=1,N
  620             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                  CALL FCN(N,X,YNEW,DY)
                  DO 621 I=1,N
  621             AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK3,IP)
                  DO 630 I=1,N
  630             YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)  
                  CALL FCN(N,X,YNEW,DY)
                  DO 631 I=1,N
  631             AK4(I)=DY(I)+HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 
              ELSE
C --- THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS
                  HD1=H*D1
                  HD2=H*D2
                  HD3=H*D3
                  HD4=H*D4
                  DO 703 I=1,N
  703             AK1(I)=DY1(I)+HD1*FX(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK1,IP)
                  DO 710 I=1,N
  710             YNEW(I)=Y(I)+A21*AK1(I) 
                  CALL FCN(N,X+C2*H,YNEW,DY)
                  DO 711 I=1,N
  711             AK2(I)=DY(I)+HD2*FX(I)+HC21*AK1(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK2,IP)
                  DO 720 I=1,N
  720             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                  CALL FCN(N,X+C3*H,YNEW,DY)
                  DO 721 I=1,N
  721             AK3(I)=DY(I)+HD3*FX(I)+HC31*AK1(I)+HC32*AK2(I)
                      CALL SOLB(N,LDE,E,MLE,MUE,AK3,IP)
                  DO 730 I=1,N
  730             YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)  
                  CALL FCN(N,X+C4*H,YNEW,DY)
                  DO 731 I=1,N
  731             AK4(I)=DY(I)+HD4*FX(I)+HC41*AK1(I)+HC42*AK2(I)
     &                  +HC43*AK3(I) 
              END IF
C --- THE DIFFERENTIAL EQUATION IS AUTONOMOUS OR NON-AUTONOMOUS
                  CALL SOLB(N,LDE,E,MLE,MUE,AK4,IP)
              DO 740 I=1,N
  740         YNEW(I)=Y(I)+A51*AK1(I)+A52*AK2(I)+A53*AK3(I)
     &                +A54*AK4(I)  
              CALL FCN(N,X+H,YNEW,DY)
              DO 741 I=1,N
  741         AK5(I)=DY(I)+HC51*AK1(I)+HC52*AK2(I)+HC53*AK3(I) 
     &              +HC54*AK4(I) 
                  CALL SOLB(N,LDE,E,MLE,MUE,AK5,IP)
C ------------ EMBEDDED SOLUTION ---------------
              DO 750 I=1,N
  750         YNEW(I)=YNEW(I)+AK5(I)  
              CALL FCN(N,X+H,YNEW,DY)
              DO 751 I=1,N
  751         AK6(I)=DY(I)+HC61*AK1(I)+HC62*AK2(I)+HC63*AK3(I)
     &              +HC64*AK4(I)+HC65*AK5(I) 
                  CALL SOLB(N,LDE,E,MLE,MUE,AK6,IP)
C ------------ NEW SOLUTION ---------------
              DO 752 I=1,N
  752         YNEW(I)=YNEW(I)+AK6(I)  
          ELSE
C + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
C --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
C ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
C ---   2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX
C + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
C --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX)
              DO 801 J=1,N
              DO 800 I=1,N
  800         E(I,J)=-FJAC(I,J)
  801         E(J,J)=E(J,J)+FAC
              CALL DEC(N,LDE,E,IP,INFO)
              IF (INFO.NE.0) GOTO 80
              IF (AUTNMS) THEN
C --- THE DIFFERENTIAL EQUATION IS AUTONOMOUS
                  DO 803 I=1,N
  803             AK1(I)=DY1(I)
                      CALL SOL(N,LDE,E,AK1,IP)
                  DO 810 I=1,N
  810             YNEW(I)=Y(I)+A21*AK1(I) 
                  CALL FCN(N,X,YNEW,DY)
                  DO 811 I=1,N
  811             AK2(I)=DY(I)+HC21*AK1(I)
                      CALL SOL(N,LDE,E,AK2,IP)
                  DO 820 I=1,N
  820             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                  CALL FCN(N,X,YNEW,DY)
                  DO 821 I=1,N
  821             AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I)
                      CALL SOL(N,LDE,E,AK3,IP)
                  DO 830 I=1,N
  830             YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)  
                  CALL FCN(N,X,YNEW,DY)
                  DO 831 I=1,N
  831             AK4(I)=DY(I)+HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) 
              ELSE
C --- THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS
                  HD1=H*D1
                  HD2=H*D2
                  HD3=H*D3
                  HD4=H*D4
                  DO 903 I=1,N
  903             AK1(I)=DY1(I)+HD1*FX(I)
                      CALL SOL(N,LDE,E,AK1,IP)
                  DO 910 I=1,N
  910             YNEW(I)=Y(I)+A21*AK1(I) 
                  CALL FCN(N,X+C2*H,YNEW,DY)
                  DO 911 I=1,N
  911             AK2(I)=DY(I)+HD2*FX(I)+HC21*AK1(I)
                      CALL SOL(N,LDE,E,AK2,IP)
                  DO 920 I=1,N
  920             YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)  
                  CALL FCN(N,X+C3*H,YNEW,DY)
                  DO 921 I=1,N
  921             AK3(I)=DY(I)+HD3*FX(I)+HC31*AK1(I)+HC32*AK2(I)
                      CALL SOL(N,LDE,E,AK3,IP)
                  DO 930 I=1,N
  930             YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I)  
                  CALL FCN(N,X+C4*H,YNEW,DY)
                  DO 931 I=1,N
  931             AK4(I)=DY(I)+HD4*FX(I)+HC41*AK1(I)+HC42*AK2(I)
     &                  +HC43*AK3(I) 
              END IF
C --- THE DIFFERENTIAL EQUATION IS AUTONOMOUS OR NON-AUTONOMOUS
                  CALL SOL(N,LDE,E,AK4,IP)
              DO 940 I=1,N
  940         YNEW(I)=Y(I)+A51*AK1(I)+A52*AK2(I)+A53*AK3(I)
     &              +A54*AK4(I)  
              CALL FCN(N,X+H,YNEW,DY)
              DO 941 I=1,N
  941         AK5(I)=DY(I)+HC51*AK1(I)+HC52*AK2(I)+HC53*AK3(I)
     &              +HC54*AK4(I) 
                  CALL SOL(N,LDE,E,AK5,IP)
C ------------ EMBEDDED SOLUTION ---------------
              DO 950 I=1,N
  950         YNEW(I)=YNEW(I)+AK5(I)  
              CALL FCN(N,X+H,YNEW,DY)
              DO 951 I=1,N
  951         AK6(I)=DY(I)+HC61*AK1(I)+HC62*AK2(I)+HC63*AK3(I)
     &              +HC64*AK4(I)+HC65*AK5(I) 
                  CALL SOL(N,LDE,E,AK6,IP)
C ------------ NEW SOLUTION ---------------
              DO 952 I=1,N
  952         YNEW(I)=YNEW(I)+AK6(I)  
          END IF
      END IF
      NSOL=NSOL+6
      NFCN=NFCN+5 
C ------------ DENSE OUTPUT ----------
      IF (IOUT.NE.0) THEN
          DO 63 I=1,N 
          CONT(I)=Y(I)
          CONT(I+NN2)=D21*AK1(I)+D22*AK2(I)+D23*AK3(I)+D24*AK4(I)
     &               +D25*AK5(I)
          CONT(I+NN3)=D31*AK1(I)+D32*AK2(I)+D33*AK3(I)+D34*AK4(I)
     &               +D35*AK5(I)
  63      CONTINUE 
      END IF
C *** *** *** *** *** *** ***
C  ERROR ESTIMATION  
C *** *** *** *** *** *** ***
      NSTEP=NSTEP+1
C ------------ COMPUTE ERROR ESTIMATION ----------------
      ERR=0.D0
      DO 300 I=1,N
      IF (ITOL.EQ.0) THEN
         SK=ATOL(1)+RTOL(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
      ELSE
         SK=ATOL(I)+RTOL(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
      END IF
  300 ERR=ERR+(AK6(I)/SK)**2
      ERR=DSQRT(ERR/N)
C --- COMPUTATION OF HNEW
C --- WE REQUIRE .2<=HNEW/H<=6.
      FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.25D0/0.9D0))
      HNEW=H/FAC  
C *** *** *** *** *** *** ***
C  IS THE ERROR SMALL ENOUGH ?
C *** *** *** *** *** *** ***
      IF (ERR.LE.1.D0) THEN
C --- STEP IS ACCEPTED  
         NACCPT=NACCPT+1
         DO 44 I=1,N 
  44     Y(I)=YNEW(I) 
         XOLD=X 
         X=X+H
         IF (IOUT.NE.0) THEN 
            DO 73 I=1,N 
  73        CONT(NN+I)=Y(I)
            IRTRN=1
            HOUT=H
            CALL SOLOUT(NACCPT+1,XOLD,X,Y,N,IRTRN)
            IF (IRTRN.LT.0) GOTO 79
         END IF
         IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN
         IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H)) 
         REJECT=.FALSE.
         H=HNEW
         GOTO 1
      ELSE
C --- STEP IS REJECTED  
         REJECT=.TRUE.
         H=HNEW
         IF (NACCPT.GE.1) NREJCT=NREJCT+1
         GOTO 2
      END IF
C --- EXIT
  80  WRITE (6,*) ' MATRIX E IS SINGULAR, INFO = ',INFO
      NSING=NSING+1
      IF (NSING.GE.5) GOTO 79
      H=H*0.5D0
      REJECT=.TRUE.
      GOTO 2
  79  WRITE(6,979)X,H
 979  FORMAT(' EXIT OF RODAS AT X=',D16.7,'   H=',D16.7)
      IDID=-1
      RETURN
      END 
C
      FUNCTION CONTRO(I,X)
C ----------------------------------------------------------
C     THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT IN CONECTION
C     WITH THE OUTPUT-SUBROUTINE FOR DOPRI5. IT PROVIDES AN
C     APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT X.
C ----------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /CONROS/N,N2,N3,XOLD,H,CONT(1)
      S=(X-XOLD)/H 
      CONTRO=CONT(I)*(1-S)+S*(CONT(I+N)+(1-S)*(CONT(I+N2)+S*CONT(I+N3)))
      RETURN
      END
C
      SUBROUTINE ROCOE1(A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,
     &  C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,
     &  C62,C63,C64,C65,GAMMA,C2,C3,C4,D1,D2,D3,D4,
     &  D21,D22,D23,D24,D25,D31,D32,D33,D34,D35)
      IMPLICIT REAL*8 (A-H,O-Z)
        C2=0.386D0
        C3=0.21D0 
        C4=0.63D0
        BET2P=0.0317D0
        BET3P=0.0635D0
        BET4P=0.3438D0 
       D1= 0.2500000000000000D+00
       D2=-0.1043000000000000D+00
       D3= 0.1035000000000000D+00
       D4=-0.3620000000000023D-01
       A21= 0.1544000000000000D+01
       A31= 0.9466785280815826D+00
       A32= 0.2557011698983284D+00
       A41= 0.3314825187068521D+01
       A42= 0.2896124015972201D+01
       A43= 0.9986419139977817D+00
       A51= 0.1221224509226641D+01
       A52= 0.6019134481288629D+01
       A53= 0.1253708332932087D+02
       A54=-0.6878860361058950D+00
       C21=-0.5668800000000000D+01
       C31=-0.2430093356833875D+01
       C32=-0.2063599157091915D+00
       C41=-0.1073529058151375D+00
       C42=-0.9594562251023355D+01
       C43=-0.2047028614809616D+02
       C51= 0.7496443313967647D+01
       C52=-0.1024680431464352D+02
       C53=-0.3399990352819905D+02
       C54= 0.1170890893206160D+02
       C61= 0.8083246795921522D+01
       C62=-0.7981132988064893D+01
       C63=-0.3152159432874371D+02
       C64= 0.1631930543123136D+02
       C65=-0.6058818238834054D+01
       GAMMA= 0.2500000000000000D+00  

       D21= 0.1012623508344586D+02
       D22=-0.7487995877610167D+01
       D23=-0.3480091861555747D+02
       D24=-0.7992771707568823D+01
       D25= 0.1025137723295662D+01
       D31=-0.6762803392801253D+00
       D32= 0.6087714651680015D+01
       D33= 0.1643084320892478D+02
       D34= 0.2476722511418386D+02
       D35=-0.6594389125716872D+01
      RETURN
      END
C
      SUBROUTINE ROCOE2(A21,A31,A32,A41,A42,A43,A51,A52,A53,A54,
     &  C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61,
     &  C62,C63,C64,C65,GAMMA,C2,C3,C4,D1,D2,D3,D4,
     &  D21,D22,D23,D24,D25,D31,D32,D33,D34,D35)
      IMPLICIT REAL*8 (A-H,O-Z)
        C2=0.3507221D0
        C3=0.2557041D0 
        C4=0.6817790D0
        BET2P=0.0317D0
        BET3P=0.0047369D0
        BET4P=0.3438D0 
       D1= 0.2500000000000000D+00
       D2=-0.6902209999999998D-01
       D3=-0.9671999999999459D-03
       D4=-0.8797900000000025D-01
       A21= 0.1402888400000000D+01
       A31= 0.6581212688557198D+00
       A32=-0.1320936088384301D+01
       A41= 0.7131197445744498D+01
       A42= 0.1602964143958207D+02
       A43=-0.5561572550509766D+01
       A51= 0.2273885722420363D+02
       A52= 0.6738147284535289D+02
       A53=-0.3121877493038560D+02
       A54= 0.7285641833203814D+00
       C21=-0.5104353600000000D+01
       C31=-0.2899967805418783D+01
       C32= 0.4040399359702244D+01
       C41=-0.3264449927841361D+02
       C42=-0.9935311008728094D+02
       C43= 0.4999119122405989D+02
       C51=-0.7646023087151691D+02
       C52=-0.2785942120829058D+03
       C53= 0.1539294840910643D+03
       C54= 0.1097101866258358D+02
       C61=-0.7629701586804983D+02
       C62=-0.2942795630511232D+03
       C63= 0.1620029695867566D+03
       C64= 0.2365166903095270D+02
       C65=-0.7652977706771382D+01
       GAMMA= 0.2500000000000000D+00  

       D21=-0.3871940424117216D+02
       D22=-0.1358025833007622D+03
       D23= 0.6451068857505875D+02
       D24=-0.4192663174613162D+01
       D25=-0.2531932050335060D+01
       D31=-0.1499268484949843D+02
       D32=-0.7630242396627033D+02
       D33= 0.5865928432851416D+02
       D34= 0.1661359034616402D+02
       D35=-0.6758691794084156D+00

      RETURN
      END
