ERG current in repolarizing plateau potentials in dopamine neurons (Canavier et al 2007)

 Download zip file 
Help downloading and running models
Accession:100603
"Blocking the small-conductance (SK) calcium-activated potassium channel promotes burst firing in dopamine neurons both in vivo and in vitro. ... We focus on the underlying plateau potential oscillation generated in the presence of both apamin and TTX, so that action potentials are not considered. We find that although the plateau potentials are mediated by a voltage-gated Ca2+ current, they do not depend on the accumulation of cytosolic Ca2+, then use a computational model to test the hypothesis that the slowly voltage-activated ether-a-go-go–related gene (ERG) potassium current repolarizes the plateaus. The model, which includes a material balance on calcium, is able to reproduce the time course of both membrane potential and somatic calcium concentration, and can also mimic the induction of plateau potentials by the calcium chelator BAPTA." See paper for more.
Reference:
1 . Canavier CC, Oprisan SA, Callaway JC, Ji H, Shepard PD (2007) Computational model predicts a role for ERG current in repolarizing plateau potentials in dopamine neurons: implications for modulation of neuronal activity. J Neurophysiol 98:3006-22 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Substantia nigra pars compacta DA cell;
Channel(s): I L high threshold; I h; I K,Ca; I Calcium; I_HERG;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program; FORTRAN;
Model Concept(s): Ion Channel Kinetics; Oscillations; Calcium dynamics;
Implementer(s): Canavier, CC;
Search NeuronDB for information about:  Substantia nigra pars compacta DA cell; I L high threshold; I h; I K,Ca; I Calcium; I_HERG;
      SUBROUTINE RADAU5(N,FCN,X,Y,XEND,H,
     &                  RTOL,ATOL,ITOL,
     &                  JAC ,IJAC,MLJAC,MUJAC,
     &                  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
C                     M*Y'=F(X,Y).
C     THE SYSTEM CAN BE (LINEARLY) IMPLICIT (MASS-MATRIX M .NE. I)
C     OR EXPLICIT (M=I).
C     THE METHOD USED IS AN IMPLICIT RUNGE-KUTTA METHOD (RADAU IIA)
C     OF ORDER 5 WITH STEP SIZE CONTROL AND CONTINUOUS OUTPUT.
C     C.F. SECTION IV.8
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 NOVEMBER 14, 1989
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     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-3 OR 1.D-5, 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 COLUMN-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)
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, RADAU5 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 REAL*8 FUNCTION
C                        >>>   CONTR5(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        SWITCH FOR CALLING THE SUBROUTINE SOLOUT:
C                    IOUT=0: SUBROUTINE IS NEVER CALLED
C                    IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT.
C
C     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
C                 WORK(1), WORK(2),.., WORK(7) SERVE AS PARAMETERS
C                 FOR THE CODE. FOR STANDARD USE OF THE CODE
C                 WORK(1),..,WORK(7) MUST BE SET TO ZERO BEFORE
C                 CALLING. SEE BELOW FOR A MORE SOPHISTICATED USE.
C                 WORK(8),..,WORK(LWORK) SERVE AS WORKING SPACE
C                 FOR ALL VECTORS AND MATRICES.
C                 "LWORK" MUST BE AT LEAST
C                             N*(LJAC+LMAS+3*LE+8)+7
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                    LE=N               IF MLJAC=N (FULL JACOBIAN)
C                    LE=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.)
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 = 4*N*N+8*N+7.
C
C     LWORK       DECLARED LENGHT OF ARRAY "WORK".
C
C     IWORK       INTEGER WORKING SPACE OF LENGHT "LIWORK".
C                 IWORK(1),IWORK(2),...,IWORK(7) SERVE AS PARAMETERS
C                 FOR THE CODE. FOR STANDARD USE, SET IWORK(1),..,
C                 IWORK(7) TO ZERO BEFORE CALLING.
C                 IWORK(8),...,IWORK(LIWORK) SERVE AS WORKING AREA.
C                 "LIWORK" MUST BE AT LEAST 3*N+7.
C
C     LIWORK      DECLARED LENGHT OF ARRAY "IWORK".
C
C     LRCONT      DECLARED LENGTH OF COMMON BLOCK
C                  >>>  COMMON /CONT/ICONT(3),RCONT(LRCONT)  <<<
C                 WHICH MUST BE DECLARED IN THE CALLING PROGRAM.
C                 "LRCONT" MUST BE AT LEAST
C                             4*N+4 .
C                 THIS IS USED FOR STORING THE COEFFICIENTS OF THE
C                 CONTINUOUS SOLUTION AND MAKES THE CALLING LIST FOR THE
C                 FUNCTION "CONTR5" 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(7)
C              AS WELL AS IWORK(1),..,IWORK(7) 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). IT IS
C              ALSO NOT GOOD FOR SPARSE JACOBIANS.
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 NEWTON ITERATIONS FOR THE
C              SOLUTION OF THE IMPLICIT SYSTEM IN EACH STEP.
C              THE DEFAULT VALUE (FOR IWORK(3)=0) IS 7.
C
C    IWORK(4)  IF IWORK(4).EQ.0 THE EXTRAPOLATED COLLOCATION SOLUTION
C              IS TAKEN AS STARTING VALUE FOR NEWTON'S METHOD.
C              IF IWORK(4).NE.0 ZERO STARTING VALUES ARE USED.
C              THE LATTER IS RECOMMENDED IF NEWTON'S METHOD HAS
C              DIFFICULTIES WITH CONVERGENCE (THIS IS THE CASE WHEN
C              NSTEP IS LARGER THAN NACCPT + NREJCT).
C              DEFAULT IS IWORK(4)=0.
C
C       THE FOLLOWING 3 PARAMETERS ARE IMPORTANT FOR
C       DIFFERENTIAL-ALGEBRAIC SYSTEMS OF INDEX > 1.
C       THE FUNCTION-SUBROUTINE SHOULD BE WRITTEN SUCH THAT
C       THE INDEX 1,2,3 VARIABLES APPEAR IN THIS ORDER.
C       IN ESTIMATING THE ERROR THE INDEX 2 VARIABLES ARE
C       MULTIPLIED BY H, THE INDEX 3 VARIABLES BY H**2.
C
C    IWORK(5)  DIMENSION OF THE INDEX 1 VARIABLES (MUST BE > 0). FOR
C              ODE'S THIS EQUALS THE DIMENSION OF THE SYSTEM.
C              DEFAULT IWORK(5)=N.
C
C    IWORK(6)  DIMENSION OF THE INDEX 2 VARIABLES. DEFAULT IWORK(6)=0.
C
C    IWORK(7)  DIMENSION OF THE INDEX 3 VARIABLES. DEFAULT IWORK(7)=0.
C
C
C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
C
C    WORK(2)   THE SAFETY FACTOR IN STEP SIZE PREDICTION,
C              DEFAULT 0.9D0.
C
C    WORK(3)   DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
C              INCREASE WORK(3), TO 0.1 SAY, WHEN JACOBIAN EVALUATIONS
C              ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER
C              (0.001D0, SAY). NEGATIV WORK(3) FORCES THE CODE THE
C              COMPUTE THE JACOBIAN AFTER EVERY ACCEPTED STEP.
C              DEFAULT 0.001D0.
C
C    WORK(4)   STOPPING CRIERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1.
C              SMALLER VALUES OF WORK(4) MAKE THE CODE SLOWER, BUT SAFER.
C              DEFAULT 0.03D0.
C
C    WORK(5) AND WORK(6) : IF WORK(5) < HNEW/HOLD < WORK(6), THEN THE
C              STEP SIZE IS NOT CHANGED. THIS SAVES, TOGETHER WITH A
C              LARGE WORK(3), LU-DECOMPOSITIONS AND COMPUTING TIME FOR
C              LARGE SYSTEMS. FOR SMALL SYSTEMS ONE MAY HAVE
C              WORK(5)=1.D0, WORK(6)=1.2D0, FOR LARGE FULL SYSTEMS
C              WORK(5)=0.99D0, WORK(6)=2.D0 MIGHT BE GOOD.
C              DEFAULTS WORK(5)=1.D0, WORK(6)=1.2D0 .
C
C    WORK(7)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
C
C-----------------------------------------------------------------------
C
C     OUTPUT PARAMETERS
C     -----------------
C     X           X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED
C                 (AFTER SUCCESSFUL RETURN X=XEND).
C
C     Y(N)        NUMERICAL SOLUTION AT X
C
C     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 IMPLCT,JBAND,ARRET,STARTN
      EXTERNAL FCN,JAC,MAS,SOLOUT
      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
C --- COMMON STAT CAN BE INSPECTED FOR STATISTICAL PURPOSES:
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 (DUE TO ERROR TEST),
C                  (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED)
C ---    NDEC      NUMBER OF LU-DECOMPOSITIONS OF BOTH MATRICES
C ---    NSOL      NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS, OF BOTH
C                  SYSTEMS; THE NSTEP FORWARD-BACKWARD SUBSTITUTIONS,
C                  NEEDED FOR STEP SIZE SELECTION, ARE NOT COUNTED
C *** *** *** *** *** *** ***
C        SETTING THE PARAMETERS
C *** *** *** *** *** *** ***
       NFCN=0
       NJAC=0
       NSTEP=0
       NACCPT=0
       NREJCT=0
       NDEC=0
       NSOL=0
       ARRET=.FALSE.
C -------- SWITCH FOR TRANSFORMATION OF JACOBIAN TO HESSIAN FORM ---
      NHESS=IWORK(1)
      IF (N.LE.2) NHESS=0
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(4,*)' WRONG INPUT IWORK(2)=',IWORK(2)
            ARRET=.TRUE.
         END IF
      END IF
C -------- NIT    MAXIMAL NUMBER OF NEWTON ITERATIONS
      IF(IWORK(3).EQ.0)THEN
         NIT=7
      ELSE
         NIT=IWORK(3)
         IF(NIT.LE.0)THEN
            WRITE(4,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
            ARRET=.TRUE.
         END IF
      END IF
C -------- STARTN  SWITCH FOR STARTING VALUES OF NEWTON ITERATIONS
      IF(IWORK(4).EQ.0)THEN
         STARTN=.FALSE.
      ELSE
         STARTN=.TRUE.
      END IF
C -------- PARAMETER FOR DIFFERENTIAL-ALGEBRAIC COMPONENTS
      NIND1=IWORK(5)
      NIND2=IWORK(6)
      NIND3=IWORK(7)
      IF (NIND1.EQ.0) NIND1=N
      IF (NIND1+NIND2+NIND3.NE.N)THEN
        WRITE(4,*)' CURIOUS INPUT FOR WORK(5,6,7)=',NIND1,NIND2,NIND3
        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.1.D-19.OR.UROUND.GE.1.D0)THEN
            WRITE(4,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',WORK(1)
            ARRET=.TRUE.
         END IF
      END IF
C --------- SAFE     SAFETY FACTOR IN STEP SIZE PREDICTION
      IF(WORK(2).EQ.0.D0)THEN
         SAFE=0.9D0
      ELSE
         SAFE=WORK(2)
         IF(SAFE.LE..001D0.OR.SAFE.GE.1.D0)THEN
            WRITE(4,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2)
            ARRET=.TRUE.
         END IF
      END IF
C ------ THET     DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
      IF(WORK(3).EQ.0.D0)THEN
         THET=0.001D0
      ELSE
         THET=WORK(3)
      END IF
C --- FNEWT   STOPPING CRIERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1.
      IF(WORK(4).EQ.0.D0)THEN
         FNEWT=0.03D0
      ELSE
         FNEWT=WORK(4)
      END IF
C --- QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST.
      IF(WORK(5).EQ.0.D0)THEN
         QUOT1=1.D0
      ELSE
         QUOT1=WORK(5)
      END IF
      IF(WORK(6).EQ.0.D0)THEN
         QUOT2=1.2D0
      ELSE
         QUOT2=WORK(6)
      END IF
C -------- MAXIMAL STEP SIZE
      IF(WORK(7).EQ.0.D0)THEN
         HMAX=XEND-X
      ELSE
         HMAX=WORK(7)
      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 ---- IMPLICIT, BANDED OR NOT ?
      IMPLCT=IMAS.NE.0
      JBAND=MLJAC.NE.N
C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
C -- JACOBIAN  AND  MATRICES E1, E2
      IF(JBAND)THEN
         LDJAC=MLJAC+MUJAC+1
         LDE1=2*MLJAC+MUJAC+1
      ELSE
         LDJAC=N
         LDE1=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 ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN
      IF ((IMPLCT.OR.JBAND).AND.NHESS.NE.0) THEN
         WRITE(4,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS WITH
     &FULL JACOBIAN'
         ARRET=.TRUE.
      END IF
C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
      IEZ1=8
      IEZ2=IEZ1+N
      IEZ3=IEZ2+N
      IEY0=IEZ3+N
      IESCAL=IEY0+N
      IEF1=IESCAL+N
      IEF2=IEF1+N
      IEF3=IEF2+N
      IEJAC=IEF3+N
      IEMAS=IEJAC+N*LDJAC
      IEE1=IEMAS+N*LDMAS
      IEE2R=IEE1+N*LDE1
      IEE2I=IEE2R+N*LDE1
C ------ TOTAL STORAGE REQUIREMENT -----------
      ISTORE=IEE2I+N*LDE1-1
      IF(ISTORE.GT.LWORK)THEN
         WRITE(4,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
         ARRET=.TRUE.
      END IF
C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
      IEIP1=8
      IEIP2=IEIP1+N
      IEIPH=IEIP2+N
C --------- TOTAL REQUIREMENT ---------------
      ISTORE=IEIPH+N-1
      IF(ISTORE.GT.LIWORK)THEN
         WRITE(4,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
         ARRET=.TRUE.
      END IF
C --------- CONTROL OF LENGTH OF COMMON BLOCK "CONT" -------
      IF(LRCONT.LT.(4*N+4))THEN
         WRITE(4,*)' INSUFF. STORAGE FOR RCONT, MIN. LRCONT=',4*N+4
         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 RADCOR(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,
     &   JAC,IJAC,MLJAC,MUJAC,MAS,IMAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
     &   NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,NHESS,STARTN,
     &   NIND1,NIND2,NIND3,
     &   IMPLCT,JBAND,LDJAC,LDE1,LDMAS2,WORK(IEZ1),WORK(IEZ2),
     &   WORK(IEZ3),WORK(IEY0),WORK(IESCAL),WORK(IEF1),WORK(IEF2),
     &   WORK(IEF3),WORK(IEJAC),WORK(IEE1),WORK(IEE2R),WORK(IEE2I),
     &   WORK(IEMAS),IWORK(IEIP1),IWORK(IEIP2),IWORK(IEIPH))
C ----------- RETURN -----------
      RETURN
      END
C
C
C
C  ----- ... AND HERE IS THE CORE INTEGRATOR  ----------
C
      SUBROUTINE RADCOR(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,
     &   JAC,IJAC,MLJAC,MUJAC,MAS,IMAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
     &   NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,NHESS,STARTN,
     &   NIND1,NIND2,NIND3,IMPLCT,BANDED,LDJAC,LDE1,LDMAS,Z1,Z2,Z3,
     &   Y0,SCAL,F1,F2,F3,FJAC,E1,E2R,E2I,FMAS,IP1,IP2,IPHES)
C ----------------------------------------------------------
C     CORE INTEGRATOR FOR RADAU5
C     PARAMETERS SAME AS IN RADAU5 WITH WORKSPACE ADDED
C ----------------------------------------------------------
C         DECLARATIONS
C ----------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 Y(N),Z1(N),Z2(N),Z3(N),Y0(N),SCAL(N),F1(N),F2(N),F3(N)
      REAL*8 FJAC(LDJAC,N),FMAS(LDMAS,N)
      REAL*8 E1(LDE1,N),E2R(LDE1,N),E2I(LDE1,N)
      REAL*8 ATOL(1),RTOL(1)
      INTEGER IP1(N),IP2(N),IPHES(N)
      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
      COMMON /CONT/NN,NN2,NN3,XSOL,HSOL,C2M1,C1M1,CONT(1)
      LOGICAL REJECT,FIRST,IMPLCT,BANDED,CALJAC,STARTN
      LOGICAL INDEX1,INDEX2,INDEX3
C *** *** *** *** *** *** ***
C  INITIALISATIONS
C *** *** *** *** *** *** ***
C --------- DUPLIFY N FOR COMMON BLOCK CONT -----
      NN=N
      NN2=2*N
      NN3=3*N
C -------- CHECK THE INDEX OF THE PROBLEM -----
      INDEX1=NIND1.NE.0
      INDEX2=NIND2.NE.0
      INDEX3=NIND3.NE.0
C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
      IF(IMPLCT)CALL MAS(N,FMAS,LDMAS)
C ---------- CONSTANTS ---------
      SQ6=DSQRT(6.D0)
      C1=(4.D0-SQ6)/10.D0
      C2=(4.D0+SQ6)/10.D0
      C1M1=C1-1.D0
      C2M1=C2-1.D0
      C1MC2=C1-C2
      DD1=-(13.D0+7.D0*SQ6)/3.D0
      DD2=(-13.D0+7.D0*SQ6)/3.D0
      DD3=-1.D0/3.D0
      U1=(6.D0+81.D0**(1.D0/3.D0)-9.D0**(1.D0/3.D0))/30.D0
      ALPH=(12.D0-81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))/60.D0
      BETA=(81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))*DSQRT(3.D0)/60.D0
      CNO=ALPH**2+BETA**2
      T11=9.1232394870892942792D-02
      T12=-0.14125529502095420843D0
      T13=-3.0029194105147424492D-02
      T21=0.24171793270710701896D0
      T22=0.20412935229379993199D0
      T23=0.38294211275726193779D0
      T31=0.96604818261509293619D0
      TI11=4.3255798900631553510D0
      TI12=0.33919925181580986954D0
      TI13=0.54177053993587487119D0
      TI21=-4.1787185915519047273D0
      TI22=-0.32768282076106238708D0
      TI23=0.47662355450055045196D0
      TI31=-0.50287263494578687595D0
      TI32=2.5719269498556054292D0
      TI33=-0.59603920482822492497D0
      POSNEG=SIGN(1.D0,XEND-X)
      HMAXN=MIN(ABS(HMAX),ABS(XEND-X))
      IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6
      H=MIN(ABS(H),HMAXN)
      H=SIGN(H,POSNEG)
      HOLD=H
      REJECT=.FALSE.
      FIRST=.TRUE.
      FACCON=1.D0
      CFAC=SAFE*(1+2*NIT)
      NSING=0
      XOLD=X
      IF (IOUT.NE.0) THEN
          IRTRN=1
          NRSOL=1
          XOSOL=XOLD
          XSOL=X
          DO 7 I=1,N
  7       CONT(I)=Y(I)
          NSOLU=N
          HSOL=HOLD
          CALL SOLOUT(NRSOL,XOSOL,XSOL,CONT,NSOLU,IRTRN)
          IF (IRTRN.LT.0) GOTO 79
      END IF
      MLE=MLJAC
      MUE=MUJAC
      MBJAC=MLJAC+MUJAC+1
      MBB=MLMAS+MUMAS+1
      MDIAG=MLE+MUE+1
      MDIFF=MLE+MUE-MUMAS
      MBDIAG=MUMAS+1
      N2=2*N
      N3=3*N
      IF (ITOL.EQ.0) THEN
          DO 8 I=1,N
   8      SCAL(I)=ATOL(1)+RTOL(1)*DABS(Y(I))
      ELSE
          DO 9 I=1,N
   9      SCAL(I)=ATOL(I)+RTOL(I)*DABS(Y(I))
      END IF
      HHFAC=H
      CALL FCN(N,X,Y,Y0)
      NFCN=NFCN+1
C --- BASIC INTEGRATION STEP
  10  CONTINUE
C *** *** *** *** *** *** ***
C  COMPUTATION OF THE JACOBIAN
C *** *** *** *** *** *** ***
      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           F1(J)=Y(J)
              F2(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J))))
              Y(J)=Y(J)+F2(J)
              J=J+MD
              IF (J.LE.N) GOTO 12
              CALL FCN(N,X,Y,CONT)
              J=K
              LBEG=MAX(1,J-MUJAC)
 14           LEND=MIN(N,J+MLJAC)
              Y(J)=F1(J)
              MUJACJ=MUJACP-J
              DO 15 L=LBEG,LEND
 15           FJAC(L+MUJACJ,J)=(CONT(L)-Y0(L))/F2(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,CONT)
              DO 17 J=1,N
  17          FJAC(J,I)=(CONT(J)-Y0(J))/DELT
  18          Y(I)=YSAFE
          END IF
      ELSE
C --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
          CALL JAC(N,X,Y,FJAC,LDJAC)
      END IF
      CALJAC=.TRUE.
      IF (NHESS.NE.0) CALL ELMHES (LDJAC,N,1,N,FJAC,IPHES)
  20  CONTINUE
C *** *** *** *** *** *** ***
C  COMPUTE THE MATRICES E1 AND E2 AND THEIR DECOMPOSITIONS
C *** *** *** *** *** *** ***
      FAC1=1.D0/(H*U1)
      ALPHN=ALPH/(H*CNO)
      BETAN=BETA/(H*CNO)
      IF (IMPLCT) THEN
          IF (BANDED) THEN
C --- THE MATRIX E (B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX)
              DO 127 J=1,N
              I1J=MAX0(1,MUJAC+2-J)
              I2J=MIN(MBJAC,MUJAC+1-J+N)
              DO 125 I=I1J,I2J
              FJ=-FJAC(I,J)
              IMLE=I+MLE
              E1(IMLE,J)=FJ
              E2R(IMLE,J)=FJ
 125          E2I(IMLE,J)=0.D0
              I1B=MAX0(1,MUMAS+2-J)
              I2B=MIN0(MBB,MUMAS+1-J+N)
              DO 126 I=I1B,I2B
              IB=I+MDIFF
              BB=FMAS(I,J)
              E1(IB,J)=E1(IB,J)+FAC1*BB
              E2R(IB,J)=E2R(IB,J)+ALPHN*BB
 126          E2I(IB,J)=BETAN*BB
 127          CONTINUE
              CALL DECB(N,LDE1,E1,MLE,MUE,IP1,IER)
              IF (IER.NE.0) GOTO 78
              CALL DECBC(N,LDE1,E2R,E2I,MLE,MUE,IP2,IER)
              IF (IER.NE.0) GOTO 78
          ELSE
              IF (MLMAS.NE.N) THEN
C --- THE MATRIX E (B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX)
                  DO 225 J=1,N
                  DO 225 I=1,N
                  FJ=-FJAC(I,J)
                  E1(I,J)=FJ
                  E2R(I,J)=FJ
 225              E2I(I,J)=0.D0
                  DO 226 J=1,N
                  I1=MAX0(1,J-MUMAS)
                  I2=MIN0(N,J+MLMAS)
                  DO 226 I=I1,I2
                  BB=FMAS(I-J+MBDIAG,J)
                  E1(I,J)=E1(I,J)+FAC1*BB
                  E2R(I,J)=E2R(I,J)+ALPHN*BB
 226              E2I(I,J)=BETAN*BB
                  CALL DEC(N,LDE1,E1,IP1,IER)
                  IF (IER.NE.0) GOTO 78
                  CALL DECC(N,LDE1,E2R,E2I,IP2,IER)
                  IF (IER.NE.0) GOTO 78
              ELSE
C --- THE MATRIX E (B IS A FULL MATRIX, JACOBIAN A FULL MATRIX)
                  IF (MLJAC.EQ.N) THEN
                      DO 324 J=1,N
                      DO 324 I=1,N
                      FJ=-FJAC(I,J)
                      BB=FMAS(I,J)
                      E1(I,J)=BB*FAC1+FJ
                      E2R(I,J)=BB*ALPHN+FJ
 324                  E2I(I,J)=BB*BETAN
                      CALL DEC(N,LDE1,E1,IP1,IER)
                      IF (IER.NE.0) GOTO 78
                      CALL DECC(N,LDE1,E2R,E2I,IP2,IER)
                      IF (IER.NE.0) GOTO 78
                  ELSE
C --- THE MATRIX E (B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX)
                  END IF
              END IF
          END IF
      ELSE
          IF (BANDED) THEN
C --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX)
              DO 427 J=1,N
              I1J=MAX0(1,MUJAC+2-J)
              I2J=MIN(MBJAC,MUJAC+1-J+N)
              DO 425 I=I1J,I2J
              FJ=-FJAC(I,J)
              IMLE=I+MLE
              E1(IMLE,J)=FJ
              E2R(IMLE,J)=FJ
 425          E2I(IMLE,J)=0.D0
              E1(MDIAG,J)=E1(MDIAG,J)+FAC1
              E2R(MDIAG,J)=E2R(MDIAG,J)+ALPHN
              E2I(MDIAG,J)=BETAN
 427          CONTINUE
              CALL DECB(N,LDE1,E1,MLE,MUE,IP1,IER)
              IF (IER.NE.0) GOTO 78
              CALL DECBC(N,LDE1,E2R,E2I,MLE,MUE,IP2,IER)
              IF (IER.NE.0) GOTO 78
          ELSE
C --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX)
              IF (NHESS.EQ.0) THEN
                  DO 526 J=1,N
                  DO 525 I=1,N
                  FJ=-FJAC(I,J)
                  E1(I,J)=FJ
                  E2R(I,J)=FJ
 525              E2I(I,J)=0.D0
                  E1(J,J)=E1(J,J)+FAC1
                  E2R(J,J)=E2R(J,J)+ALPHN
 526              E2I(J,J)=BETAN
                  CALL DEC(N,LDE1,E1,IP1,IER)
                  IF (IER.NE.0) GOTO 78
                  CALL DECC(N,LDE1,E2R,E2I,IP2,IER)
                  IF (IER.NE.0) GOTO 78
              ELSE
                  DO 624 J=1,N-1
                  J1=J+1
                  FJ=-FJAC(J1,J)
                  E1(J1,J)=FJ
                  E2R(J1,J)=FJ
 624              E2I(J1,J)=0.D0
                  DO 626 J=1,N
                  DO 625 I=1,J
                  FJ=-FJAC(I,J)
                  E1(I,J)=FJ
                  E2I(I,J)=0.D0
 625              E2R(I,J)=FJ
                  E1(J,J)=E1(J,J)+FAC1
                  E2R(J,J)=E2R(J,J)+ALPHN
 626              E2I(J,J)=BETAN
                  CALL DECH(N,LDE1,E1,1,IP1,IER)
                  IF (IER.NE.0) GOTO 78
                  CALL DECHC(N,LDE1,E2R,E2I,1,IP2,IER)
                  IF (IER.NE.0) GOTO 78
              END IF
          END IF
      END IF
      NDEC=NDEC+1
  30  CONTINUE
      NSTEP=NSTEP+1
      IF (NSTEP.GT.NMAX.OR.X+1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79
          IF (INDEX2) THEN
             DO 465 I=NIND1+1,NIND1+NIND2
 465         SCAL(I)=SCAL(I)/HHFAC
          END IF
          IF (INDEX3) THEN
             DO 466 I=NIND1+NIND2+1,NIND1+NIND2+NIND3
 466         SCAL(I)=SCAL(I)/(HHFAC*HHFAC)
          END IF
      XPH=X+H
C *** *** *** *** *** *** ***
C  STARTING VALUES FOR NEWTON ITERATION
C *** *** *** *** *** *** ***
      IF (FIRST.OR.STARTN) THEN
         DO 32 I=1,N
         Z1(I)=0.D0
         Z2(I)=0.D0
         Z3(I)=0.D0
         F1(I)=0.D0
         F2(I)=0.D0
  32     F3(I)=0.D0
      ELSE
         C3Q=H/HOLD
         C1Q=C1*C3Q
         C2Q=C2*C3Q
         DO 35 I=1,N
         AK1=CONT(I+N)
         AK2=CONT(I+N2)
         AK3=CONT(I+N3)
         Z1I=C1Q*(AK1+(C1Q-C2M1)*(AK2+(C1Q-C1M1)*AK3))
         Z2I=C2Q*(AK1+(C2Q-C2M1)*(AK2+(C2Q-C1M1)*AK3))
         Z3I=C3Q*(AK1+(C3Q-C2M1)*(AK2+(C3Q-C1M1)*AK3))
         Z1(I)=Z1I
         Z2(I)=Z2I
         Z3(I)=Z3I
         F1(I)=TI11*Z1I+TI12*Z2I+TI13*Z3I
         F2(I)=TI21*Z1I+TI22*Z2I+TI23*Z3I
  35     F3(I)=TI31*Z1I+TI32*Z2I+TI33*Z3I
      END IF
C *** *** *** *** *** *** ***
C  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
C *** *** *** *** *** *** ***
            NEWT=0
            FACCON=MAX(FACCON,UROUND)**0.8D0
            THETA=ABS(THET)
  40        CONTINUE
            IF (NEWT.GE.NIT) GOTO 78
C ---     COMPUTE THE RIGHT-HAND SIDE
            DO 41 I=1,N
  41        CONT(I)=Y(I)+Z1(I)
            CALL FCN(N,X+C1*H,CONT,Z1)
            DO 42 I=1,N
  42        CONT(I)=Y(I)+Z2(I)
            CALL FCN(N,X+C2*H,CONT,Z2)
            DO 43 I=1,N
  43        CONT(I)=Y(I)+Z3(I)
            CALL FCN(N,XPH,CONT,Z3)
            NFCN=NFCN+3
C ---     SOLVE THE LINEAR SYSTEMS
            IF (IMPLCT) THEN
                IF (MLMAS.NE.N) THEN
                    DO 146 I=1,N
                    S1=0.0D0
                    S2=0.0D0
                    S3=0.0D0
                    J1B=MAX0(1,I-MLMAS)
                    J2B=MIN0(N,I+MUMAS)
                    DO 145 J=J1B,J2B
                    BB=FMAS(I-J+MBDIAG,J)
                    S1=S1-BB*F1(J)
                    S2=S2-BB*F2(J)
 145                S3=S3-BB*F3(J)
                    A1=Z1(I)
                    A2=Z2(I)
                    A3=Z3(I)
                    Z1(I)=S1*FAC1          +TI11*A1+TI12*A2+TI13*A3
                    Z2(I)=S2*ALPHN-S3*BETAN+TI21*A1+TI22*A2+TI23*A3
 146                Z3(I)=S3*ALPHN+S2*BETAN+TI31*A1+TI32*A2+TI33*A3
                    IF (BANDED) THEN
                        CALL SOLB(N,LDE1,E1,MLE,MUE,Z1,IP1)
                        CALL SOLBC(N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
                    ELSE
                        CALL SOL(N,LDE1,E1,Z1,IP1)
                        CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
                    END IF
                ELSE
                    DO 246 I=1,N
                    S1=0.0D0
                    S2=0.0D0
                    S3=0.0D0
                    DO 245 J=1,N
                    BB=FMAS(I,J)
                    S1=S1-BB*F1(J)
                    S2=S2-BB*F2(J)
 245                S3=S3-BB*F3(J)
                    A1=Z1(I)
                    A2=Z2(I)
                    A3=Z3(I)
                    Z1(I)=S1*FAC1          +TI11*A1+TI12*A2+TI13*A3
                    Z2(I)=S2*ALPHN-S3*BETAN+TI21*A1+TI22*A2+TI23*A3
 246                Z3(I)=S3*ALPHN+S2*BETAN+TI31*A1+TI32*A2+TI33*A3
                    CALL SOL(N,LDE1,E1,Z1,IP1)
                    CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
                END IF
            ELSE
                DO 345 I=1,N
                A1=Z1(I)
                A2=Z2(I)
                A3=Z3(I)
                S2=-F2(I)
                S3=-F3(I)
                Z1(I)=-F1(I)*FAC1      +TI11*A1+TI12*A2+TI13*A3
                Z2(I)=S2*ALPHN-S3*BETAN+TI21*A1+TI22*A2+TI23*A3
 345            Z3(I)=S3*ALPHN+S2*BETAN+TI31*A1+TI32*A2+TI33*A3
                IF (BANDED) THEN
                    CALL SOLB(N,LDE1,E1,MLE,MUE,Z1,IP1)
                    CALL SOLBC(N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
                ELSE
                    IF (NHESS.EQ.0) THEN
                        CALL SOL(N,LDE1,E1,Z1,IP1)
                        CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
                    ELSE
                        DO 140 MM=N-2,1,-1
                        MP=N-MM
                        MP1=MP-1
                        I=IPHES(MP)
                        IF (I.EQ.MP) GOTO 110
                        ZSAFE=Z1(MP)
                        Z1(MP)=Z1(I)
                        Z1(I)=ZSAFE
                        ZSAFE=Z2(MP)
                        Z2(MP)=Z2(I)
                        Z2(I)=ZSAFE
                        ZSAFE=Z3(MP)
                        Z3(MP)=Z3(I)
                        Z3(I)=ZSAFE
 110                    CONTINUE
                        DO 100 I=MP+1,N
                        E1IMP=FJAC(I,MP1)
                        Z1(I)=Z1(I)-E1IMP*Z1(MP)
                        Z2(I)=Z2(I)-E1IMP*Z2(MP)
 100                    Z3(I)=Z3(I)-E1IMP*Z3(MP)
 140                    CONTINUE
                           CALL SOLH(N,LDE1,E1,1,Z1,IP1)
                           CALL SOLHC(N,LDE1,E2R,E2I,1,Z2,Z3,IP2)
                        DO 240 MM=1,N-2
                        MP=N-MM
                        MP1=MP-1
                        DO 200 I=MP+1,N
                        E1IMP=FJAC(I,MP1)
                        Z1(I)=Z1(I)+E1IMP*Z1(MP)
                        Z2(I)=Z2(I)+E1IMP*Z2(MP)
 200                    Z3(I)=Z3(I)+E1IMP*Z3(MP)
                        I=IPHES(MP)
                        IF (I.EQ.MP) GOTO 240
                        ZSAFE=Z1(MP)
                        Z1(MP)=Z1(I)
                        Z1(I)=ZSAFE
                        ZSAFE=Z2(MP)
                        Z2(MP)=Z2(I)
                        Z2(I)=ZSAFE
                        ZSAFE=Z3(MP)
                        Z3(MP)=Z3(I)
                        Z3(I)=ZSAFE
 240                    CONTINUE
                    END IF
                END IF
            END IF
            DO 51 I=1,N
            F1(I)=F1(I)+Z1(I)
            F2(I)=F2(I)+Z2(I)
  51        F3(I)=F3(I)+Z3(I)
            NSOL=NSOL+1
            NEWT=NEWT+1
            DYNO=0.D0
            DO 57 I=1,N
            DENOM=SCAL(I)
  57        DYNO=DYNO+(Z1(I)/DENOM)**2+(Z2(I)/DENOM)**2
     &          +(Z3(I)/DENOM)**2
            DYNO=DSQRT(DYNO/N3)
C ---     BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE
            IF (NEWT.GE.2.AND.NEWT.LT.NIT) THEN
                THETA=DYNO/DYNOLD
                IF (THETA.LT.0.99D0) THEN
                    FACCON=THETA/(1.0D0-THETA)
                    DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)
                    QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH/FNEWT))
                    IF (QNEWT.GE.1.0D0) THEN
                         HHFAC=.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT))
                         H=HHFAC*H
                         REJECT=.TRUE.
                         IF (CALJAC) GOTO 20
                         GOTO 10
                    END IF
                ELSE
                    GOTO 78
                END IF
            END IF
            DYNOLD=DMAX1(DYNO,UROUND)
            DO 58 I=1,N
            F1I=F1(I)
            F2I=F2(I)
            F3I=F3(I)
            Z1(I)=T11*F1I+T12*F2I+T13*F3I
            Z2(I)=T21*F1I+T22*F2I+T23*F3I
  58        Z3(I)=T31*F1I+    F2I
            IF (FACCON*DYNO.GT.FNEWT) GOTO 40
C *** *** *** *** *** *** ***
C  ERROR ESTIMATION
C *** *** *** *** *** *** ***
      HEE1=DD1/H
      HEE2=DD2/H
      HEE3=DD3/H
      IF (IMPLCT) THEN
          DO 359 I=1,N
 359      F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
          IF (MLMAS.EQ.N) THEN
              DO 361 I=1,N
              SUM=0.D0
              DO 360 J=1,N
 360          SUM=SUM+FMAS(I,J)*F1(J)
 361          F2(I)=SUM
          ELSE
              DO 363 I=1,N
              SUM=0.D0
              J1B=MAX0(1,I-MLMAS)
              J2B=MIN0(N,I+MUMAS)
              DO 362 J=J1B,J2B
 362          SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J)
 363          F2(I)=SUM
          END IF
      ELSE
          DO 461 I=1,N
 461      F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
      END IF
      DO 462 I=1,N
 462  CONT(I)=F2(I)+Y0(I)
      IF (BANDED) THEN
          CALL SOLB(N,LDE1,E1,MLE,MUE,CONT,IP1)
      ELSE
          IF (NHESS.EQ.0) THEN
              CALL SOL(N,LDE1,E1,CONT,IP1)
          ELSE
              DO 340 MM=N-2,1,-1
              MP=N-MM
              I=IPHES(MP)
              IF (I.EQ.MP) GOTO 310
              ZSAFE=CONT(MP)
              CONT(MP)=CONT(I)
              CONT(I)=ZSAFE
 310          CONTINUE
              DO 300 I=MP+1,N
 300          CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
 340          CONTINUE
                CALL SOLH(N,LDE1,E1,1,CONT,IP1)
              DO 440 MM=1,N-2
              MP=N-MM
              DO 400 I=MP+1,N
 400          CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
              I=IPHES(MP)
              IF (I.EQ.MP) GOTO 440
              ZSAFE=CONT(MP)
              CONT(MP)=CONT(I)
              CONT(I)=ZSAFE
 440          CONTINUE
          END IF
      END IF
      ERR=0.D0
         DO 354 I=1,N
 354     ERR=ERR+(CONT(I)/SCAL(I))**2
      ERR=DMAX1(DSQRT(ERR/N),1.D-10)
      IF ((FIRST.OR.REJECT).AND.ERR.GE.1.D0) THEN
          DO 460 I=1,N
 460      CONT(I)=Y(I)+CONT(I)
          CALL FCN(N,X,CONT,F1)
          NFCN=NFCN+1
          DO 463 I=1,N
 463      CONT(I)=F1(I)+F2(I)
          IF (BANDED) THEN
              CALL SOLB(N,LDE1,E1,MLE,MUE,CONT,IP1)
          ELSE
              IF (NHESS.EQ.0) THEN
                  CALL SOL(N,LDE1,E1,CONT,IP1)
              ELSE
                  DO 540 MM=N-2,1,-1
                  MP=N-MM
                  I=IPHES(MP)
                  IF (I.EQ.MP) GOTO 510
                  ZSAFE=CONT(MP)
                  CONT(MP)=CONT(I)
                  CONT(I)=ZSAFE
 510              CONTINUE
                  DO 500 I=MP+1,N
 500              CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
 540              CONTINUE
                    CALL SOLH(N,LDE1,E1,1,CONT,IP1)
                  DO 640 MM=1,N-2
                  MP=N-MM
                  DO 600 I=MP+1,N
 600              CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
                  I=IPHES(MP)
                  IF (I.EQ.MP) GOTO 640
                  ZSAFE=CONT(MP)
                  CONT(MP)=CONT(I)
                  CONT(I)=ZSAFE
 640              CONTINUE
              END IF
          END IF
          ERR=0.D0
             DO 364 I=1,N
 364         ERR=ERR+(CONT(I)/SCAL(I))**2
          ERR=DMAX1(DSQRT(ERR/N),1.D-10)
      END IF
C --- COMPUTATION OF HNEW
C --- WE REQUIRE .2<=HNEW/H<=8.
      FAC=DMIN1(SAFE,CFAC/(NEWT+2*NIT))
      QUOT=DMAX1(.125D0,DMIN1(5.D0,ERR**.25D0/FAC))
      HNEW=H/QUOT
C *** *** *** *** *** *** ***
C  IS THE ERROR SMALL ENOUGH ?
C *** *** *** *** *** *** ***
      IF (ERR.LT.1.D0) THEN
C --- STEP IS ACCEPTED
         FIRST=.FALSE.
         NACCPT=NACCPT+1
         XOLD=X
         HOLD=H
         X=XPH
         DO 75 I=1,N
         Y(I)=Y(I)+Z3(I)
         Z2I=Z2(I)
         Z1I=Z1(I)
         CONT(I+N)=(Z2I-Z3(I))/C2M1
         AK=(Z1I-Z2I)/C1MC2
         ACONT3=Z1I/C1
         ACONT3=(AK-ACONT3)/C2
         CONT(I+N2)=(AK-CONT(I+N))/C1M1
  75     CONT(I+N3)=CONT(I+N2)-ACONT3
         IF (ITOL.EQ.0) THEN
             DO 81 I=1,N
  81         SCAL(I)=ATOL(1)+RTOL(1)*DABS(Y(I))
         ELSE
             DO 82 I=1,N
  82         SCAL(I)=ATOL(I)+RTOL(I)*DABS(Y(I))
         END IF
         IF (IOUT.NE.0) THEN
             NRSOL=NACCPT+1
             XSOL=X
             XOSOL=XOLD
             DO 77 I=1,N
  77         CONT(I)=Y(I)
             NSOLU=N
             HSOL=HOLD
             CALL SOLOUT(NRSOL,XOSOL,XSOL,CONT,NSOLU,IRTRN)
             IF (IRTRN.LT.0) GOTO 79
         END IF
         CALJAC=.FALSE.
         IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN
            H=HOPT
            IDID=1
            RETURN
         END IF
         CALL FCN(N,X,Y,Y0)
         NFCN=NFCN+1
         HNEW=POSNEG*DMIN1(DABS(HNEW),HMAXN)
         HOPT=HNEW
         IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H))
         REJECT=.FALSE.
         IF ((X+HNEW/QUOT1-XEND)*POSNEG.GT.0.D0) THEN
            H=XEND-X
         ELSE
            QT=HNEW/H
            HHFAC=H
            IF (THETA.LE.THET.AND.QT.GE.QUOT1.AND.QT.LE.QUOT2) GOTO 30
            H=HNEW
         END IF
         HHFAC=H
         IF (THETA.LE.THET) GOTO 20
         GOTO 10
      ELSE
C --- STEP IS REJECTED
         REJECT=.TRUE.
         IF (FIRST) THEN
             H=H*0.1D0
             HHFAC=0.1D0
         ELSE
             HHFAC=HNEW/H
             H=HNEW
         END IF
         IF (NACCPT.GE.1) NREJCT=NREJCT+1
         IF (CALJAC) GOTO 20
         GOTO 10
      END IF
C --- UNEXPECTED STEP-REJECTION
  78  CONTINUE
      IF (IER.NE.0) THEN
          WRITE (6,*) ' MATRIX IS SINGULAR, IER=',IER,' X=',X,' H=',H
          NSING=NSING+1
          IF (NSING.GE.6) GOTO 79
      END IF
      H=H*0.5D0
      HHFAC=0.5D0
      REJECT=.TRUE.
      IF (CALJAC) GOTO 20
      GOTO 10
C --- FAIL EXIT
  79  WRITE (6,*) X,H,IER,NSTEP,IRTRN,NSING
 979  FORMAT('X=',D14.7,'   H=',D14.7,'   IER=',I4,'   NSTEP=')
      IDID=-1
      RETURN
      END
C
C
      REAL*8 FUNCTION CONTR5(I,X)
C ----------------------------------------------------------
C     THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT. IT PROVIDES AN
C     APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT X.
C     IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR
C     THE LAST SUCCESSFULLY COMPUTED STEP (BY RADAU5).
C ----------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /CONT/NN,NN2,NN3,XSOL,HSOL,C2M1,C1M1,CONT(1)
      S=(X-XSOL)/HSOL
      CONTR5=CONT(I)+S*(CONT(I+NN)+(S-C2M1)*(CONT(I+NN2)
     &     +(S-C1M1)*CONT(I+NN3)))
      RETURN
      END
C



Loading data, please wait...