SUBROUTINE SEULEX(N,FCN,IFCN,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JAC ,IJAC,MLJAC,MUJAC,
     &                  MAS ,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,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 EXTRAPOLATION-ALGORITHM, BASED ON THE
C     LINEARLY IMPLICIT EULER METHOD (WITH STEP SIZE CONTROL
C     AND ORDER SELECTION).
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              INCLUSION OF DENSE OUTPUT BY E. HAIRER AND A. OSTERMANN
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 14,
C         SPRINGER-VERLAG 1991, SECOND EDITION 1996.
C      
C     VERSION OF 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                 RPAR, IPAR (SEE BELOW)
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 (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,RPAR,IPAR)
C                    DOUBLE PRECISION 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     ----   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,RPAR,IPAR)
C                    DOUBLE PRECISION 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,RC,LRC,IC,LIC,N,
C                                       RPAR,IPAR,IRTRN)
C                    DOUBLE PRECISION X,Y(N),RC(LRC),IC(LIC)
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, SEULEX RETURNS TO THE CALLING PROGRAM.
C                 DO NOT CHANGE THE ENTRIES OF RC(LRC),IC(LIC)!
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                        >>>   CONTEX(I,S,RC,LRC,IC,LIC)   <<<
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                    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 AND MATRICES.
C                 "LWORK" MUST BE AT LEAST
C                        N*(LJAC+LMAS+LE1+KM+8)+4*KM+20+KM2*NRDENS
C                 WHERE
C                    KM2=2+KM*(KM+3)/2  AND  NRDENS=IWORK(6) (SEE BELOW)
C                 AND
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                 AND                                 
C                    KM=12               IF IWORK(3)=0
C                    KM=IWORK(3)         IF IWORK(3).GT.0
C
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+(KM+8)*N+4*KM+13+KM2*NRDENS.
C                 IF IWORK(9)=M1>0 THEN "LWORK" MUST BE AT LEAST
C                    N*(LJAC+KM+8)+(N-M1)*(LMAS+LE1)+4*KM+20+KM2*NRDENS
C                 WHERE IN THE DEFINITIONS OF LJAC, LMAS AND LE1 THE
C                 NUMBER N CAN BE REPLACED BY N-M1.
C
C     LWORK       DECLARED LENGTH OF ARRAY "WORK".
C
C     IWORK       INTEGER WORKING SPACE OF LENGTH "LIWORK".
C                 "LIWORK" MUST BE AT LEAST  2*N+KM+20+NRDENS.
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 OF THE CODE ARE TUNED TO MAKE IT WORK 
C              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(13)
C              AS WELL AS IWORK(1),..,IWORK(4) DIFFERENT FROM ZERO.
C              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
C
C    IWORK(1)  IF IWORK(1).NE.0, THE CODE TRANSFORMS THE JACOBIAN
C              MATRIX TO HESSENBERG FORM. THIS IS PARTICULARLY
C              ADVANTAGEOUS FOR LARGE SYSTEMS WITH FULL JACOBIAN.
C              IT DOES NOT WORK FOR BANDED JACOBIAN (MLJAC<N)
C              AND NOT FOR IMPLICIT SYSTEMS (IMAS=1).
C
C    IWORK(2)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
C              THE DEFAULT VALUE (FOR IWORK(2)=0) IS 100000.
C
C    IWORK(3)  THE MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION 
C              TABLE. THE DEFAULT VALUE (FOR IWORK(3)=0) IS 12.
C              IF IWORK(3).NE.0 THEN IWORK(3) SHOULD BE .GE.3.
C
C    IWORK(4)  SWITCH FOR THE STEP SIZE SEQUENCE
C              IF IWORK(4).EQ.1 THEN 1,2,3,4,6,8,12,16,24,32,48,...
C              IF IWORK(4).EQ.2 THEN 2,3,4,6,8,12,16,24,32,48,64,...
C              IF IWORK(4).EQ.3 THEN 1,2,3,4,5,6,7,8,9,10,...
C              IF IWORK(4).EQ.4 THEN 2,3,4,5,6,7,8,9,10,11,... 
C              THE DEFAULT VALUE (FOR IWORK(4)=0) IS IWORK(4)=2.
C
C    IWORK(5)  PARAMETER "LAMBDA" OF DENSE OUTPUT; POSSIBLE VALUES
C              ARE 0 AND 1; DEFAULT IWORK(5)=0.
C
C    IWORK(6)  = 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
C
C       IF THE DIFFERENTIAL SYSTEM HAS THE SPECIAL STRUCTURE THAT
C            Y(I)' = Y(I+M2)   FOR  I=1,...,M1,
C       WITH M1 A MULTIPLE OF M2, A SUBSTANTIAL GAIN IN COMPUTERTIME
C       CAN BE ACHIEVED BY SETTING THE PARAMETERS IWORK(9) AND IWORK(10).
C       E.G., FOR SECOND ORDER SYSTEMS P'=V, V'=G(P,V), WHERE P AND V ARE 
C       VECTORS OF DIMENSION N/2, ONE HAS TO PUT M1=M2=N/2.
C       FOR M1>0 SOME OF THE INPUT PARAMETERS HAVE DIFFERENT MEANINGS:
C       - JAC: ONLY THE ELEMENTS OF THE NON-TRIVIAL PART OF THE
C              JACOBIAN HAVE TO BE STORED
C              IF (MLJAC.EQ.N-M1) THE JACOBIAN IS SUPPOSED TO BE FULL
C                 DFY(I,J) = PARTIAL F(I+M1) / PARTIAL Y(J)
C                FOR I=1,N-M1 AND J=1,N.
C              ELSE, THE JACOBIAN IS BANDED ( M1 = M2 * MM )
C                 DFY(I-J+MUJAC+1,J+K*M2) = PARTIAL F(I+M1) / PARTIAL Y(J+K*M2)
C                FOR I=1,MLJAC+MUJAC+1 AND J=1,M2 AND K=0,MM.
C       - MLJAC: MLJAC=N-M1: IF THE NON-TRIVIAL PART OF THE JACOBIAN IS FULL
C                0<=MLJAC<N-M1: IF THE (MM+1) SUBMATRICES (FOR K=0,MM)
C                     PARTIAL F(I+M1) / PARTIAL Y(J+K*M2),  I,J=1,M2
C                    ARE BANDED, MLJAC IS THE MAXIMAL LOWER BANDWIDTH
C                    OF THESE MM+1 SUBMATRICES
C       - MUJAC: MAXIMAL UPPER BANDWIDTH OF THESE MM+1 SUBMATRICES
C                NEED NOT BE DEFINED IF MLJAC=N-M1
C       - MAS: IF IMAS=0 THIS MATRIX IS ASSUMED TO BE THE IDENTITY AND
C              NEED NOT BE DEFINED. SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
C              IT IS ASSUMED THAT ONLY THE ELEMENTS OF RIGHT LOWER BLOCK OF
C              DIMENSION N-M1 DIFFER FROM THAT OF THE IDENTITY MATRIX.
C              IF (MLMAS.EQ.N-M1) THIS SUBMATRIX IS SUPPOSED TO BE FULL
C                 AM(I,J) = M(I+M1,J+M1)     FOR I=1,N-M1 AND J=1,N-M1.
C              ELSE, THE MASS MATRIX IS BANDED
C                 AM(I-J+MUMAS+1,J) = M(I+M1,J+M1)
C       - MLMAS: MLMAS=N-M1: IF THE NON-TRIVIAL PART OF M IS FULL
C                0<=MLMAS<N-M1: LOWER BANDWIDTH OF THE MASS MATRIX
C       - MUMAS: UPPER BANDWIDTH OF THE MASS MATRIX
C                NEED NOT BE DEFINED IF MLMAS=N-M1
C
C    IWORK(9)  THE VALUE OF M1.  DEFAULT M1=0.
C
C    IWORK(10) THE VALUE OF M2.  DEFAULT M2=M1.
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)   DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
C              INCREASE WORK(3), TO 0.01 SAY, WHEN JACOBIAN EVALUATIONS
C              ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER.
C              DEFAULT MIN(1.0D-4,RTOL(1))
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/(J-1)) 
C              DEFAULT VALUES: WORK(4)=0.1D0, WORK(5)=4.D0
C
C    WORK(6), WORK(7)   PARAMETERS FOR THE ORDER SELECTION
C              ORDER IS DECREASED IF    W(K-1) <= W(K)*WORK(6)
C              ORDER IS INCREASED IF    W(K) <= W(K-1)*WORK(7)
C              DEFAULT VALUES: WORK(6)=0.7D0, 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.8D0, WORK(9)=0.93D0
C
C    WORK(10), WORK(11), WORK(12), WORK(13)   ESTIMATED WORKS FOR
C             A CALL TO  FCN, JAC, DEC, SOL, RESPECTIVELY.
C             DEFAULT VALUES ARE: WORK(10)=1.D0, WORK(11)=5.D0,
C             WORK(12)=1.D0, WORK(13)=1.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   IWORK(14)  NFCN    NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL
C                      EVALUATION OF THE JACOBIAN ARE NOT COUNTED)  
C   IWORK(15)  NJAC    NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
C                      OR NUMERICALLY)
C   IWORK(16)  NSTEP   NUMBER OF COMPUTED STEPS
C   IWORK(17)  NACCPT  NUMBER OF ACCEPTED STEPS
C   IWORK(18)  NREJCT  NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP
C                      HAS BEEN ACCEPTED)
C   IWORK(19)  NDEC    NUMBER OF LU-DECOMPOSITIONS (N-DIMENSIONAL MATRIX)
C   IWORK(20)  NSOL    NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS
C-----------------------------------------------------------------------
C *** *** *** *** *** *** *** *** *** *** *** *** ***
C          DECLARATIONS 
C *** *** *** *** *** *** *** *** *** *** *** *** ***
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(N),ATOL(*),RTOL(*),WORK(LWORK),IWORK(LIWORK)
      DIMENSION RPAR(*),IPAR(*)
      LOGICAL AUTNMS,IMPLCT,ARRET,JBAND
      EXTERNAL FCN,JAC,MAS,SOLOUT
C *** *** *** *** *** *** ***
C        SETTING THE PARAMETERS 
C *** *** *** *** *** *** ***
      NFCN=0
      NJAC=0
      NSTEP=0
      NACCPT=0
      NREJCT=0
      NDEC=0
      NSOL=0
      ARRET=.FALSE.
C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
      IF(IWORK(2).EQ.0)THEN
         NMAX=100000
      ELSE
         NMAX=IWORK(2)
         IF(NMAX.LE.0)THEN
            WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2)
            ARRET=.TRUE.
         END IF
      END IF
C -------- KM     MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION 
      IF(IWORK(3).EQ.0)THEN
         KM=12
      ELSE
         KM=IWORK(3)
         IF(KM.LE.2)THEN
            WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
            ARRET=.TRUE.
         END IF
      END IF
C -------- NSEQU     CHOICE OF STEP SIZE SEQUENCE
      NSEQU=IWORK(4)
      IF(IWORK(4).EQ.0) NSEQU=2
      IF(NSEQU.LE.0.OR.NSEQU.GE.5)THEN
         WRITE(6,*)' CURIOUS INPUT IWORK(4)=',IWORK(4)
         ARRET=.TRUE.
      END IF 
C -------- LAMBDA   PARAMETER FOR DENSE OUTPUT
      LAMBDA=IWORK(5)
      IF(LAMBDA.LT.0.OR.LAMBDA.GE.2)THEN
         WRITE(6,*)' CURIOUS INPUT IWORK(5)=',IWORK(5)
         ARRET=.TRUE.
      END IF 
C -------- NRDENS   NUMBER OF DENSE OUTPUT COMPONENTS
      NRDENS=IWORK(6)
      IF(NRDENS.LT.0.OR.NRDENS.GT.N)THEN
         WRITE(6,*)' CURIOUS INPUT IWORK(6)=',IWORK(6)
         ARRET=.TRUE.
      END IF
C -------- PARAMETER FOR SECOND ORDER EQUATIONS
      M1=IWORK(9)
      M2=IWORK(10)
      NM1=N-M1
      IF (M1.EQ.0) M2=N
      IF (M2.EQ.0) M2=M1
      IF (M1.LT.0.OR.M2.LT.0.OR.M1+M2.GT.N) THEN
       WRITE(6,*)' CURIOUS INPUT FOR IWORK(9,10)=',M1,M2
       ARRET=.TRUE.
      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.0.D0.OR.UROUND.GE.1.D0)THEN
            WRITE(6,*)'  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 ------ THET     DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
      IF(WORK(3).EQ.0.D0)THEN
         THET=MIN(1.0D-4,RTOL(1))
      ELSE
         THET=WORK(3)
      END IF
C -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
      IF(WORK(4).EQ.0.D0)THEN
         FAC1=0.1D0
      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.7D0
      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.6D0
      ELSE
         SAFE1=WORK(8)
      END IF
      IF(WORK(9).EQ.0.D0)THEN
         SAFE2=0.93D0
      ELSE
         SAFE2=WORK(9)
      END IF
C ------- WKFCN,WKJAC,WKDEC,WKSOL  ESTIMATED WORK FOR  FCN,JAC,DEC,SOL
      IF(WORK(10).EQ.0.D0)THEN
         WKFCN=1.D0
      ELSE
         WKFCN=WORK(10)
      END IF
      IF(WORK(11).EQ.0.D0)THEN
         WKJAC=5.D0
      ELSE
         WKJAC=WORK(11)
      END IF
      IF(WORK(12).EQ.0.D0)THEN
         WKDEC=1.D0
      ELSE
         WKDEC=WORK(12)
      END IF
      IF(WORK(13).EQ.0.D0)THEN
         WKSOL=1.D0
      ELSE
         WKSOL=WORK(13)
      END IF
      WKROW=WKFCN+WKSOL
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 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
        END DO
      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.LT.NM1
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
         MLJAC=NM1
         MUJAC=NM1
         LDJAC=NM1
         LDE=NM1
      END IF
C -- MASS MATRIX
      IF (IMPLCT) THEN
          IF (MLMAS.NE.NM1) THEN
              LDMAS=MLMAS+MUMAS+1
              IF (JBAND) THEN
                 IJOB=4
              ELSE
                 IJOB=3
              END IF
          ELSE
              LDMAS=NM1
              IJOB=5
          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
          IF (JBAND) THEN
             IJOB=2
          ELSE
             IJOB=1
             IF (N.GT.2.AND.IWORK(1).NE.0) IJOB=7
          END IF
      END IF
      LDMAS2=MAX(1,LDMAS)
C ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN
      IF ((IMPLCT.OR.JBAND).AND.IJOB.EQ.7) THEN
         WRITE(6,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS WITH 
     &FULL JACOBIAN'
         ARRET=.TRUE.
      END IF
C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- 
      KM2=(KM*(KM+1))/2
      IEYH=21
      IEDY=IEYH+N
      IEFX=IEDY+N
      IEYHH=IEFX+N
      IEDYH=IEYHH+N
      IEDEL=IEDYH+N
      IEWH =IEDEL+N
      IESCAL=IEWH+N
      IEHH =IESCAL+N
      IEW  =IEHH+KM
      IEA  =IEW+KM
      IEJAC =IEA+KM
      IEE  =IEJAC+N*LDJAC
      IEMAS=IEE+NM1*LDE
      IET=IEMAS+NM1*LDMAS
      IFAC=IET+N*KM  
      IEDE=IFAC+KM
      IFSAFE=IEDE+(KM+2)*NRDENS
C ------ TOTAL STORAGE REQUIREMENT -----------
      ISTORE=IFSAFE+KM2*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 -----
      IECO=21
      IEIP=21+NRDENS
      IENJ=IEIP+N 
      IEIPH=IENJ+KM
C --------- TOTAL REQUIREMENT ---------------
      ISTORE=IECO+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 
      NRD=MAX(1,NRDENS)
C -------- CALL TO CORE INTEGRATOR ------------
      CALL SEUCOR(N,FCN,X,Y,XEND,HMAX,H,KM,RTOL,ATOL,ITOL,JAC,IJAC,
     &   MLJAC,MUJAC,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,IJOB,M1,M2,NM1,
     &   NMAX,UROUND,NSEQU,AUTNMS,IMPLCT,JBAND,LDJAC,LDE,LDMAS2,
     &   WORK(IEYH),WORK(IEDY),WORK(IEFX),WORK(IEYHH),WORK(IEDYH),
     &   WORK(IEDEL),WORK(IEWH),WORK(IESCAL),WORK(IEHH),
     &   WORK(IEW),WORK(IEA),WORK(IEJAC),WORK(IEE),WORK(IEMAS),
     &   WORK(IET),IWORK(IEIP),IWORK(IENJ),IWORK(IEIPH),FAC1,FAC2,FAC3,
     &   FAC4,THET,SAFE1,SAFE2,WKJAC,WKDEC,WKROW,KM2,NRD,WORK(IFAC),
     &   WORK(IFSAFE),LAMBDA,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,
     &   RPAR,IPAR,WORK(IEDE),IWORK(IECO))
      IWORK(14)=NFCN
      IWORK(15)=NJAC
      IWORK(16)=NSTEP
      IWORK(17)=NACCPT
      IWORK(18)=NREJCT
      IWORK(19)=NDEC
      IWORK(20)=NSOL
C ----------- RETURN -----------
      RETURN
      END
C
C
C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
C
      SUBROUTINE SEUCOR(N,FCN,X,Y,XEND,HMAX,H,KM,RTOL,ATOL,ITOL,
     &  JAC,IJAC,MLJAC,MUJAC,MAS,MLB,MUB,SOLOUT,IOUT,IDID,IJOB,M1,M2,
     &  NM1,NMAX,UROUND,NSEQU,AUTNMS,IMPLCT,BANDED,LFJAC,LE,
     &  LDMAS,YH,DY,FX,YHH,DYH,DEL,WH,SCAL,HH,W,A,FJAC,E,FMAS,T,IP,
     &  NJ,IPHES,FAC1,FAC2,FAC3,FAC4,THET,SAFE1,SAFE2,WKJAC,WKDEC,WKROW,
     &  KM2,NRD,FACUL,FSAFE,LAMBDA,NFCN,NJAC,NSTEP,NACCPT,NREJCT,
     &  NDEC,NSOL,RPAR,IPAR,DENS,ICOMP)
C ----------------------------------------------------------
C     CORE INTEGRATOR FOR SEULEX
C     PARAMETERS SAME AS IN SEULEX WITH WORKSPACE ADDED 
C ---------------------------------------------------------- 
C         DECLARATIONS 
C ---------------------------------------------------------- 
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DIMENSION Y(N),YH(N),DY(N),FX(N),YHH(N),DYH(N),DEL(N)
       DIMENSION WH(N),SCAL(N),HH(KM),W(KM),A(KM),FJAC(LFJAC,N)
       DIMENSION FMAS(LDMAS,NM1),T(KM,N),IP(NM1),NJ(KM),RTOL(*),ATOL(*)
       DIMENSION IPHES(N),FSAFE(KM2,NRD),FACUL(KM),E(LE,NM1)
       DIMENSION DENS((KM+2)*NRD),ICOMP(NRD),RPAR(*),IPAR(*)
       LOGICAL REJECT,LAST,ATOV,CALJAC,CALHES,AUTNMS,IMPLCT,BANDED
       EXTERNAL FCN
       COMMON /COSEU/XOLDD,HHH,NNRD,KRIGHT
       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C --- COMPUTE COEFFICIENTS FOR DENSE OUTPUT
       IF (IOUT.EQ.2) THEN
          NNRD=NRD
C --- COMPUTE THE FACTORIALS --------
          FACUL(1)=1.D0
          DO I=1,KM-1
             FACUL(I+1)=I*FACUL(I)
          END DO
       END IF
C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
      IF (IMPLCT) CALL MAS(NM1,FMAS,LDMAS,RPAR,IPAR)
C *** *** *** *** *** *** ***
C  INITIALISATIONS
C *** *** *** *** *** *** ***
       LRDE=(KM+2)*NRD
       MLE=MLJAC
       MUE=MUJAC
       MBJAC=MLJAC+MUJAC+1
       MBB=MLB+MUB+1
       MDIAG=MLE+MUE+1
       MDIFF=MLE+MUE-MUB
       MBDIAG=MUB+1
       IF (M1.GT.0) IJOB=IJOB+10
C --- DEFINE THE STEP SIZE SEQUENCE
       IF (NSEQU.EQ.1) THEN
           NJ(1)=1
           NJ(2)=2
           NJ(3)=3
           DO I=4,KM
              NJ(I)=2*NJ(I-2)
           END DO
       END IF
       IF (NSEQU.EQ.2) THEN
           NJ(1)=2
           NJ(2)=3
           DO I=3,KM
              NJ(I)=2*NJ(I-2)
           END DO
       END IF
       DO I=1,KM
          IF (NSEQU.EQ.3) NJ(I)=I
          IF (NSEQU.EQ.4) NJ(I)=I+1
       END DO
       A(1)=WKJAC+NJ(1)*WKROW+WKDEC
       DO I=2,KM
          A(I)=A(I-1)+(NJ(I)-1)*WKROW+WKDEC
       END DO
       K=MAX0(3,MIN0(KM-2,INT(-DLOG10(RTOL(1)+ATOL(1))*.6D0+1.5D0))) 
       HMAXN=MIN(ABS(HMAX),ABS(XEND-X)) 
       IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6
       H=MIN(ABS(H),HMAXN)
       THETA=2*ABS(THET)
       IF (IOUT.NE.0) THEN
          IRTRN=1
          NRSOL=1
          XSOL=X
          DO I=1,N
            YH(I)=Y(I)
          END DO
          NSOLU=N
          XOLD=X
         CALL SOLOUT(NRSOL,XOLD,XSOL,YH,DENS,LRDE,ICOMP,NRD,NSOLU,
     &               RPAR,IPAR,IRTRN)
          IF (IRTRN.LT.0) GOTO 120
          IPT=0
       END IF
       ERR=0.D0
       W(1)=1.D30  
       DO I=1,N
          IF (ITOL.EQ.0) THEN
            SCAL(I)=ATOL(1)+RTOL(1)*DABS(Y(I))
          ELSE
            SCAL(I)=ATOL(I)+RTOL(I)*DABS(Y(I))
          END IF
       END DO
       CALJAC=.FALSE.
       REJECT=.FALSE.
       LAST=.FALSE.
  10   CONTINUE
       IF (REJECT) THETA=2*ABS(THET)
       ATOV=.FALSE.
C *** *** *** *** *** *** ***
C --- IS XEND REACHED IN THE NEXT STEP?
C *** *** *** *** *** *** ***
       H1=XEND-X  
       IF (H1.LE.UROUND) GO TO 110  
       HOPT=H
       H=MIN(H,H1,HMAXN)
       IF (H.GE.H1-UROUND) LAST=.TRUE.
       IF (AUTNMS) THEN   
          CALL FCN(N,X,Y,DY,RPAR,IPAR)
          NFCN=NFCN+1
       END IF
       IF (THETA.GT.THET.AND..NOT.CALJAC) THEN
C *** *** *** *** *** *** ***
C  COMPUTATION OF THE JACOBIAN
C *** *** *** *** *** *** ***
          NJAC=NJAC+1
          IF (IJAC.EQ.0) THEN
C --- COMPUTE JACOBIAN MATRIX NUMERICALLY
              IF (.NOT.AUTNMS) CALL FCN(N,X,Y,DY,RPAR,IPAR)
              IF (BANDED) THEN
C --- JACOBIAN IS BANDED
                  MUJACP=MUJAC+1
                  MD=MIN(MBJAC,N)
                  DO MM=1,M1/M2+1
                     DO KX=1,MD
                        J=KX+(MM-1)*M2
 12                     YHH(J)=Y(J)
                        DEL(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J))))
                        Y(J)=Y(J)+DEL(J)
                        J=J+MD
                        IF (J.LE.MM*M2) GOTO 12 
                        CALL FCN(N,X,Y,YH,RPAR,IPAR)
                        J=KX+(MM-1)*M2
                        J1=KX
                        LBEG=MAX(1,J1-MUJAC)+M1
 14                     LEND=MIN(M2,J1+MLJAC)+M1
                        Y(J)=YHH(J)
                        MUJACJ=MUJACP-J1-M1
                        DO L=LBEG,LEND
                           FJAC(L+MUJACJ,J)=(YH(L)-DY(L))/DEL(J) 
                        END DO
                        J=J+MD
                        J1=J1+MD
                        LBEG=LEND+1
                        IF (J.LE.MM*M2) GOTO 14
                     END DO
                  END DO
               ELSE
C --- JACOBIAN IS FULL
                  DO I=1,N
                     YSAFE=Y(I)
                     DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE)))
                     Y(I)=YSAFE+DELT
                     CALL FCN(N,X,Y,YH,RPAR,IPAR)
                     DO J=M1+1,N
                       FJAC(J-M1,I)=(YH(J)-DY(J))/DELT
                     END DO
                     Y(I)=YSAFE
                     MLJAC=NM1
                  END DO
               END IF
            ELSE
C --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
               CALL JAC(N,X,Y,FJAC,LFJAC,RPAR,IPAR)
           END IF
           CALJAC=.TRUE.
           CALHES=.TRUE.
       END IF
C *** *** *** *** *** *** ***
C --- THE FIRST AND LAST STEP 
C *** *** *** *** *** *** ***
      IF (NSTEP.EQ.0.OR.LAST) THEN 
        IPT=0
        NSTEP=NSTEP+1 
        DO J=1,K
         KC=J
         CALL SEUL(J,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM,
     &             HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1,FAC,
     &             FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL,MLB,
     &             MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT,
     &             ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES,
     &             RPAR,IPAR)
         IF (ATOV) GOTO 10
         IF (J.GT.1.AND.ERR.LE.1.D0) GO TO 60
        END DO
        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 J=1,KC
       CALL SEUL(J,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM,
     &           HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1,
     &           FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL,
     &           MLB,MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT,
     &           ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES,
     &           RPAR,IPAR)
       IF (ATOV) GOTO 10
      END DO
C *** *** *** *** *** *** ***
C --- CONVERGENCE MONITOR
C *** *** *** *** *** *** ***
      IF (K.EQ.2.OR.REJECT) GO TO 50
      IF (ERR.LE.1.D0) GO TO 60
      IF (ERR.GT.DBLE(NJ(K+1)*NJ(K))*4.D0) GO TO 100  
  50  CALL SEUL(K,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM,
     &           HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1,
     &           FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL,
     &           MLB,MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT,
     &           ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES,
     &           RPAR,IPAR)
      IF (ATOV) GOTO 10
      KC=K 
      IF (ERR.LE.1.D0) GO TO 60
C --- HOPE FOR CONVERGENCE IN LINE K+1
  55  IF (ERR.GT.DBLE(NJ(K+1))*2.D0) GO TO 100  
      KC=K+1
      CALL SEUL(KC,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM,
     &           HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1,
     &           FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL,
     &           MLB,MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT,
     &           ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES,
     &           RPAR,IPAR)
      IF (ATOV) GOTO 10
      IF (ERR.GT.1.D0) GO TO 100
C *** *** *** *** *** *** ***
C --- STEP IS ACCEPTED  
C *** *** *** *** *** *** ***
  60  XOLD=X 
      X=X+H
      IF (IOUT.EQ.2) THEN 
         KRIGHT=KC
         DO I=1,NRD
            DENS(I)=Y(ICOMP(I))
         END DO
      END IF
      DO I=1,N
         T1I=T(1,I)
         IF (ITOL.EQ.0) THEN
            SCAL(I)=ATOL(1)+RTOL(1)*DABS(T1I)
         ELSE
            SCAL(I)=ATOL(I)+RTOL(I)*DABS(T1I)
         END IF
         Y(I)=T1I
      END DO
      NACCPT=NACCPT+1  
      CALJAC=.FALSE.  
      IF (IOUT.EQ.2) THEN
         XOLDD=XOLD
         HHH=H
         DO I=1,NRD
            DENS(NRD+I)=Y(ICOMP(I))
         END DO
         DO KLR=1,KRIGHT-1 
C --- COMPUTE DIFFERENCES
            IF (KLR.GE.2) THEN  
               DO KK=KLR,KC
                  LBEG=((KK+1)*KK)/2 
                  LEND=LBEG-KK+2
                  DO L=LBEG,LEND,-1
                     DO I=1,NRD
                        FSAFE(L,I)=FSAFE(L,I)-FSAFE(L-1,I)
                     END DO
                  END DO
               END DO
             END IF
C --- COMPUTE DERIVATIVES AT RIGHT END ----
             DO KK=KLR+LAMBDA,KC
                FACNJ=NJ(KK)
                FACNJ=FACNJ**KLR/FACUL(KLR+1)
                IPT=((KK+1)*KK)/2
                DO  I=1,NRD   
                   KRN=(KK-LAMBDA+1)*NRD
                   DENS(KRN+I)=FSAFE(IPT,I)*FACNJ
                END DO
             END DO
             DO J=KLR+LAMBDA+1,KC
               DBLENJ=NJ(J)
               DO L=J,KLR+LAMBDA+1,-1
                 FACTOR=DBLENJ/NJ(L-1)-1.D0
                 DO I=1,NRD
                   KRN=(L-LAMBDA+1)*NRD+I
                   DENS(KRN-NRD)=DENS(KRN)
     &                          +(DENS(KRN)-DENS(KRN-NRD))/FACTOR
                 END DO  
               END DO  
             END DO  
         END DO  
C ---  COMPUTE THE COEFFICIENTS OF THE INTERPOLATION POLYNOMIAL
         DO IN=1,NRD
            DO J=1,KRIGHT 
               II=NRD*J+IN
               DENS(II)=DENS(II)-DENS(II-NRD)
            END DO
         END DO
      END IF
      IF (IOUT.NE.0) THEN
         IRTRN=1
         NRSOL=NACCPT+1
         XSOL=X
         DO I=1,N
            YH(I)=Y(I)
         END DO
         NSOLU=N
         CALL SOLOUT(NRSOL,XOLD,XSOL,YH,DENS,LRDE,ICOMP,NRD,NSOLU,
     &               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=MIN0(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=MIN0(KC,KM-1)
      END IF
C --- AFTER A REJECTED STEP
  80  IF (REJECT) THEN 
         K=MIN0(KOPT,KC)
         H=DMIN1(H,HH(K))
         REJECT=.FALSE.
         GO TO 10
      END IF
C --- COMPUTE STEP SIZE 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
      GO TO 10
C *** *** *** *** *** *** ***
C --- STEP IS REJECTED  
C *** *** *** *** *** *** ***
 100  K=MIN0(K,KC)
      IF (K.GT.2.AND.W(K-1).LT.W(K)*FAC3) K=K-1
      NREJCT=NREJCT+1  
      H=HH(K)
      LAST=.FALSE.
      REJECT=.TRUE.
      IF (CALJAC) GOTO 30
      GO TO 10
C --- SOLUTION EXIT
 110  CONTINUE
      H=HOPT
      IDID=1
      RETURN
C --- FAIL EXIT
 120  WRITE (6,979) X,H
 979  FORMAT(' EXIT OF SEULEX AT X=',D14.7,'   H=',D14.7)
      IDID=-1
      RETURN
      END  
C
C
C *** *** *** *** *** *** ***
C     S U B R O U T I N E    S E U L
C *** *** *** *** *** *** ***
C
      SUBROUTINE SEUL(JJ,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,
     &          H,KM,HMAXN,T,SCAL,NJ,HH,W,A,YH,DYH,DEL,WH,ERR,SAFE1,
     &          FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL,
     &          MLB,MUB,ERROLD,IPHES,ICOMP,
     &          AUTNMS,IMPLCT,BANDED,REJECT,ATOV,FSAFE,KM2,NRD,IOUT,
     &          IPT,M1,M2,NM1,IJOB,CALHES,RPAR,IPAR)
C --- THIS SUBROUTINE COMPUTES THE J-TH LINE OF THE
C --- EXTRAPOLATION TABLE AND PROVIDES AN ESTIMATE  
C --- OF THE OPTIMAL STEP SIZE 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(N),YH(N),DY(N),FX(N),DYH(N),DEL(N)
      DIMENSION WH(N),SCAL(N),HH(KM),W(KM),A(KM),FJAC(LFJAC,N),E(LE,N)
      DIMENSION FMAS(LDMAS,N),T(KM,N),IP(N),NJ(KM),IPHES(N)
      DIMENSION FSAFE(KM2,NRD),ICOMP(NRD),RPAR(*),IPAR(*)
      LOGICAL ATOV,REJECT,AUTNMS,IMPLCT,BANDED,CALHES
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C *** *** *** *** *** *** ***
C  COMPUTE THE MATRIX E AND ITS DECOMPOSITION
C *** *** *** *** *** *** ***
      HJ=H/NJ(JJ)
      HJI=1.D0/HJ
      CALL DECOMR(N,FJAC,LFJAC,FMAS,LDMAS,MLB,MUB,
     &            M1,M2,NM1,HJI,E,LE,IP,IER,IJOB,CALHES,IPHES)
      IF (IER.NE.0) GOTO 79
      NDEC=NDEC+1
C *** *** *** *** *** *** ***
C --- STARTING PROCEDURE
C *** *** *** *** *** *** ***
       IF (.NOT.AUTNMS) THEN
          CALL FCN(N,X+HJ,Y,DY,RPAR,IPAR)
          NFCN=NFCN+1   
       END IF
       DO I=1,N
          YH(I)=Y(I)
          DEL(I)=DY(I)
       END DO
       CALL SLVSEU(N,FJAC,LFJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
     &          M1,M2,NM1,HJI,E,LE,IP,IPHES,DEL,IJOB)
       NSOL=NSOL+1
       M=NJ(JJ)
       IF (IOUT.EQ.2.AND.M.EQ.JJ) THEN 
          IPT=IPT+1
          DO I=1,NRD
             FSAFE(IPT,I)=DEL(ICOMP(I))
          END DO
       END IF
C *** *** *** *** *** *** ***
C --- SEMI-IMPLICIT EULER METHOD
C *** *** *** *** *** *** ***
       IF (M.GT.1) THEN
          DO MM=1,M-1
             DO I=1,N
                YH(I)=YH(I)+DEL(I)
             END DO
             IF (AUTNMS) THEN
                CALL FCN(N,X+HJ*MM,YH,DYH,RPAR,IPAR)
             ELSE
                CALL FCN(N,X+HJ*(MM+1),YH,DYH,RPAR,IPAR)
             END IF
             NFCN=NFCN+1
             IF (MM.EQ.1.AND.JJ.LE.2) THEN
C --- STABILITY CHECK
                DEL1=0.D0
                DO I=1,N
                   DEL1=DEL1+(DEL(I)/SCAL(I))**2
                END DO
                DEL1=DSQRT(DEL1)
                IF (IMPLCT) THEN
                   DO I=1,NM1
                      WH(I)=DEL(I+M1)
                   END DO
                   IF (MLB.EQ.NM1) THEN
                      DO I=1,NM1
                         SUM=0.D0
                         DO J=1,NM1
                            SUM=SUM+FMAS(I,J)*WH(J)
                         END DO
                         DEL(I+M1)=SUM
                      END DO
                   ELSE
                      DO I=1,NM1
                         SUM=0.D0
                         DO J=MAX(1,I-MLB),MIN(NM1,I+MUB)
                            SUM=SUM+FMAS(I-J+MBDIAG,J)*WH(J)
                         END DO
                         DEL(I+M1)=SUM
                      END DO
                   END IF
                END IF
                IF (.NOT.AUTNMS) THEN
                   CALL FCN(N,X+HJ,YH,WH,RPAR,IPAR)  
                   NFCN=NFCN+1
                   DO I=1,N
                      DEL(I)=WH(I)-DEL(I)*HJI
                   END DO
                ELSE
                   DO I=1,N
                      DEL(I)=DYH(I)-DEL(I)*HJI
                   END DO
                END IF
                CALL SLVSEU(N,FJAC,LFJAC,MLJAC,MUJAC,FMAS,LDMAS,
     &            MLMAS,MUMAS,M1,M2,NM1,HJI,E,LE,IP,IPHES,DEL,IJOB)
                NSOL=NSOL+1
                DEL2=0.D0
                DO I=1,N
                   DEL2=DEL2+(DEL(I)/SCAL(I))**2
                END DO
                DEL2=DSQRT(DEL2)
                THETA=DEL2/MAX(1.D0,DEL1)
                IF (THETA.GT.1.D0) GOTO 79
             END IF
             CALL SLVSEU(N,FJAC,LFJAC,MLJAC,MUJAC,FMAS,LDMAS,
     &          MLMAS,MUMAS,M1,M2,NM1,HJI,E,LE,IP,IPHES,DYH,IJOB)
             NSOL=NSOL+1
             DO I=1,N
                DEL(I)=DYH(I)
             END DO
             IF (IOUT.EQ.2.AND.MM.GE.M-JJ) THEN 
                IPT=IPT+1
                DO I=1,NRD
                   FSAFE(IPT,I)=DEL(ICOMP(I))
                END DO
             END IF
          END DO
       END IF
       DO I=1,N
          T(JJ,I)=YH(I)+DEL(I)
       END DO
C *** *** *** *** *** *** ***
C --- POLYNOMIAL EXTRAPOLATION
C *** *** *** *** *** *** ***
       IF (JJ.EQ.1) RETURN
       DO L=JJ,2,-1
          FAC=(DBLE(NJ(JJ))/DBLE(NJ(L-1)))-1.D0
          DO I=1,N
             T(L-1,I)=T(L,I)+(T(L,I)-T(L-1,I))/FAC
          END DO
       END DO
       ERR=0.D0
       DO I=1,N
          ERR=ERR+MIN(ABS((T(1,I)-T(2,I)))/SCAL(I),1.D15)**2 
       END DO  
       IF (ERR.GE.1.D30) GOTO 79
       ERR=DSQRT(ERR/DBLE(N))
       IF (JJ.GT.2.AND.ERR.GE.ERROLD) GOTO 79
       ERROLD=DMAX1(4*ERR,1.D0)
C --- COMPUTE OPTIMAL STEP SIZES
       EXPO=1.D0/JJ
       FACMIN=FAC1**EXPO
       FAC=DMIN1(FAC2/FACMIN,DMAX1(FACMIN,(ERR/SAFE1)**EXPO/SAFE2))
       FAC=1.D0/FAC
       HH(JJ)=DMIN1(H*FAC,HMAXN)
       W(JJ)=A(JJ)/HH(JJ)  
       RETURN
  79   ATOV=.TRUE.
       H=H*0.5D0
       REJECT=.TRUE.
       RETURN
       END 
C
      FUNCTION CONTEX(II,X,RC,LRC,IC,LIC)
C ----------------------------------------------------------
C     THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT IN CONECTION
C     WITH THE OUTPUT-SUBROUTINE FOR SEULEX. IT PROVIDES AN
C     APPROXIMATION TO THE II-TH COMPONENT OF THE SOLUTION AT X.
C ----------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION RC(LRC),IC(LIC)
      COMMON /COSEU/XOLD,H,NRD,IR
C ----- COMPUTE PLACE OF II-TH COMPONENT 
      I=0 
      DO J=1,NRD
        IF (IC(J).EQ.II) I=J
      END DO
      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 
      CONTEX=RC(IR*NRD+I)
      DO J=2,IR
         CONTEX=RC((1+IR-J)*NRD+I)+CONTEX*(THETA-1.D0)
      END DO
      CONTEX=RC(I)+CONTEX*THETA
      RETURN
      END