1      SUBROUTINE RADAU(N,FCN,X,Y,XEND,H,
2     &                  RTOL,ATOL,ITOL,
3     &                  JAC ,IJAC,MLJAC,MUJAC,
4     &                  MAS ,IMAS,MLMAS,MUMAS,
5     &                  SOLOUT,IOUT,
6     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
7C ----------------------------------------------------------
8C     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
9C     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS
10C                     M*Y'=F(X,Y).
11C     THE SYSTEM CAN BE (LINEARLY) IMPLICIT (MASS-MATRIX M .NE. I)
12C     OR EXPLICIT (M=I).
13C     THE CODE IS BASED ON IMPLICIT RUNGE-KUTTA METHODS (RADAU IIA)
14C     WITH VARIABLE ORDER (5, 9, 13), WITH STEP SIZE CONTROL
15C     AND CONTINUOUS OUTPUT.
16C
17C     AUTHORS: E. HAIRER AND G. WANNER
18C              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
19C              CH-1211 GENEVE 24, SWITZERLAND
20C              E-MAIL:  Ernst.Hairer@math.unige.ch
21C                       Gerhard.Wanner@math.unige.ch
22C
23C     FOR A DESCRIPTION OF THE RELATED CODE RADAU5 SEE THE BOOK:
24C         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
25C         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
26C         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14,
27C         SPRINGER-VERLAG 1991, SECOND EDITION 1996.
28C
29C     PRELIMINARY VERSION OF APRIL 23, 1998
30C     (latest small correction: January 18, 2002)
31C
32C     INPUT PARAMETERS
33C     ----------------
34C     N           DIMENSION OF THE SYSTEM
35C
36C     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
37C                 VALUE OF F(X,Y):
38C                    SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR)
39C                    DOUBLE PRECISION X,Y(N),F(N)
40C                    F(1)=...   ETC.
41C                 RPAR, IPAR (SEE BELOW)
42C
43C     X           INITIAL X-VALUE
44C
45C     Y(N)        INITIAL VALUES FOR Y
46C
47C     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
48C
49C     H           INITIAL STEP SIZE GUESS;
50C                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,
51C                 H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD.
52C                 THIS CHOICE IS NOT VERY IMPORTANT, THE STEP SIZE IS
53C                 QUICKLY ADAPTED. (IF H=0.D0, THE CODE PUTS H=1.D-6).
54C
55C     RTOL,ATOL   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
56C                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
57C
58C     ITOL        SWITCH FOR RTOL AND ATOL:
59C                   ITOL=0: BOTH RTOL AND ATOL ARE SCALARS.
60C                     THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF
61C                     Y(I) BELOW RTOL*ABS(Y(I))+ATOL
62C                   ITOL=1: BOTH RTOL AND ATOL ARE VECTORS.
63C                     THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW
64C                     RTOL(I)*ABS(Y(I))+ATOL(I).
65C
66C     JAC         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
67C                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y
68C                 (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY
69C                 A DUMMY SUBROUTINE IN THE CASE IJAC=0).
70C                 FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM
71C                    SUBROUTINE JAC(N,X,Y,DFY,LDFY,RPAR,IPAR)
72C                    DOUBLE PRECISION X,Y(N),DFY(LDFY,N)
73C                    DFY(1,1)= ...
74C                 LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS
75C                 FURNISHED BY THE CALLING PROGRAM.
76C                 IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO
77C                    BE FULL AND THE PARTIAL DERIVATIVES ARE
78C                    STORED IN DFY AS
79C                       DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J)
80C                 ELSE, THE JACOBIAN IS TAKEN AS BANDED AND
81C                    THE PARTIAL DERIVATIVES ARE STORED
82C                    DIAGONAL-WISE AS
83C                       DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).
84C
85C     IJAC        SWITCH FOR THE COMPUTATION OF THE JACOBIAN:
86C                    IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE
87C                       DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED.
88C                    IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.
89C
90C     MLJAC       SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN:
91C                    MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR
92C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
93C                    0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN
94C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
95C                       THE MAIN DIAGONAL).
96C
97C     MUJAC       UPPER BANDWITH OF JACOBIAN  MATRIX (>= NUMBER OF NON-
98C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
99C                 NEED NOT BE DEFINED IF MLJAC=N.
100C
101C     ----   MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS      -----
102C     ----   FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): -
103C
104C     MAS         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS-
105C                 MATRIX M.
106C                 IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY
107C                 MATRIX AND NEEDS NOT TO BE DEFINED;
108C                 SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
109C                 IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM
110C                    SUBROUTINE MAS(N,AM,LMAS,RPAR,IPAR)
111C                    DOUBLE PRECISION AM(LMAS,N)
112C                    AM(1,1)= ....
113C                    IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED
114C                    AS FULL MATRIX LIKE
115C                         AM(I,J) = M(I,J)
116C                    ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED
117C                    DIAGONAL-WISE AS
118C                         AM(I-J+MUMAS+1,J) = M(I,J).
119C
120C     IMAS       GIVES INFORMATION ON THE MASS-MATRIX:
121C                    IMAS=0: M IS SUPPOSED TO BE THE IDENTITY
122C                       MATRIX, MAS IS NEVER CALLED.
123C                    IMAS=1: MASS-MATRIX  IS SUPPLIED.
124C
125C     MLMAS       SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX:
126C                    MLMAS=N: THE FULL MATRIX CASE. THE LINEAR
127C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
128C                    0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE
129C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
130C                       THE MAIN DIAGONAL).
131C                 MLMAS IS SUPPOSED TO BE .LE. MLJAC.
132C
133C     MUMAS       UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON-
134C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
135C                 NEED NOT BE DEFINED IF MLMAS=N.
136C                 MUMAS IS SUPPOSED TO BE .LE. MUJAC.
137C
138C     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
139C                 NUMERICAL SOLUTION DURING INTEGRATION.
140C                 IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
141C                 SUPPLY A DUMMY SUBROUTINE IF IOUT=0.
142C                 IT MUST HAVE THE FORM
143C                    SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,
144C                                       RPAR,IPAR,IRTRN)
145C                    DOUBLE PRECISION X,Y(N),CONT(LRC)
146C                    ....
147C                 SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
148C                    GRID-POINT "X" (THEREBY THE INITIAL VALUE IS
149C                    THE FIRST GRID-POINT).
150C                 "XOLD" IS THE PRECEEDING GRID-POINT.
151C                 "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
152C                    IS SET <0, RADAU RETURNS TO THE CALLING PROGRAM.
153C
154C          -----  CONTINUOUS OUTPUT: -----
155C                 DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
156C                 FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH
157C                 THE FUNCTION
158C                        >>>   CONTRA(I,S,CONT,LRC)   <<<
159C                 WHICH PROVIDES AN APPROXIMATION TO THE I-TH
160C                 COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
161C                 S SHOULD LIE IN THE INTERVAL [XOLD,X].
162C
163C     IOUT        SWITCH FOR CALLING THE SUBROUTINE SOLOUT:
164C                    IOUT=0: SUBROUTINE IS NEVER CALLED
165C                    IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT.
166C
167C     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
168C                 WORK(1), WORK(2),.., WORK(20) SERVE AS PARAMETERS
169C                 FOR THE CODE. FOR STANDARD USE OF THE CODE
170C                 WORK(1),..,WORK(20) MUST BE SET TO ZERO BEFORE THE
171C                 FIRST CALL. SEE BELOW FOR A MORE SOPHISTICATED USE.
172C                 WORK(8),..,WORK(LWORK) SERVE AS WORKING SPACE
173C                 FOR ALL VECTORS AND MATRICES.
174C                 "LWORK" MUST BE AT LEAST
175C                          N*(LJAC+LMAS+NSMAX*LE+3*NSMAX+3)+20
176C                 WHERE
177C                    NSMAX=IWORK(12) (SEE BELOW)
178C                 AND
179C                    LJAC=N              IF MLJAC=N (FULL JACOBIAN)
180C                    LJAC=MLJAC+MUJAC+1  IF MLJAC<N (BANDED JAC.)
181C                 AND
182C                    LMAS=0              IF IMAS=0
183C                    LMAS=N              IF IMAS=1 AND MLMAS=N (FULL)
184C                    LMAS=MLMAS+MUMAS+1  IF MLMAS<N (BANDED MASS-M.)
185C                 AND
186C                    LE=N               IF MLJAC=N (FULL JACOBIAN)
187C                    LE=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.)
188C
189C                 IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE
190C                 MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM
191C                 STORAGE REQUIREMENT IS
192C                      LWORK = (NSMAX+1)*N*N+(3*NSMAX+3)*N+20.
193C                 IF IWORK(9)=M1>0 THEN "LWORK" MUST BE AT LEAST
194C                      N*(LJAC+3*NSMAX+3)+(N-M1)*(LMAS+NSMAX*LE)+20
195C                 WHERE IN THE DEFINITIONS OF LJAC, LMAS AND LE THE
196C                 NUMBER N CAN BE REPLACED BY N-M1.
197C
198C     LWORK       DECLARED LENGTH OF ARRAY "WORK".
199C
200C     IWORK       INTEGER WORKING SPACE OF LENGTH "LIWORK".
201C                 IWORK(1),IWORK(2),...,IWORK(20) SERVE AS PARAMETERS
202C                 FOR THE CODE. FOR STANDARD USE, SET IWORK(1),..,
203C                 IWORK(20) TO ZERO BEFORE CALLING.
204C                 IWORK(21),...,IWORK(LIWORK) SERVE AS WORKING AREA.
205C                 "LIWORK" MUST BE AT LEAST
206C                             (2+(NSMAX-1)/2)*N+20.
207C
208C     LIWORK      DECLARED LENGTH OF ARRAY "IWORK".
209C
210C     RPAR, IPAR  REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH
211C                 CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING
212C                 PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES.
213C
214C ----------------------------------------------------------------------
215C
216C     SOPHISTICATED SETTING OF PARAMETERS
217C     -----------------------------------
218C              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK
219C              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),...
220C              AS WELL AS IWORK(1),... DIFFERENT FROM ZERO.
221C              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
222C
223C    IWORK(1)  IF IWORK(1).NE.0, THE CODE TRANSFORMS THE JACOBIAN
224C              MATRIX TO HESSENBERG FORM. THIS IS PARTICULARLY
225C              ADVANTAGEOUS FOR LARGE SYSTEMS WITH FULL JACOBIAN.
226C              IT DOES NOT WORK FOR BANDED JACOBIAN (MLJAC<N)
227C              AND NOT FOR IMPLICIT SYSTEMS (IMAS=1).
228C
229C    IWORK(2)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
230C              THE DEFAULT VALUE (FOR IWORK(2)=0) IS 100000.
231C
232C    IWORK(3)  THE MAXIMUM NUMBER OF NEWTON ITERATIONS FOR THE
233C              SOLUTION OF THE IMPLICIT SYSTEM IN EACH STEP
234C              IWORK(3)+(NS-3)*2.5. DEFAULT VALUE (FOR IWORK(3)=0) IS 7.
235C              NS IS THE NUMBER OF STAGES (SEE IWORK(11)).
236C
237C    IWORK(4)  IF IWORK(4).EQ.0 THE EXTRAPOLATED COLLOCATION SOLUTION
238C              IS TAKEN AS STARTING VALUE FOR NEWTON'S METHOD.
239C              IF IWORK(4).NE.0 ZERO STARTING VALUES ARE USED.
240C              THE LATTER IS RECOMMENDED IF NEWTON'S METHOD HAS
241C              DIFFICULTIES WITH CONVERGENCE (THIS IS THE CASE WHEN
242C              NSTEP IS LARGER THAN NACCPT + NREJCT; SEE OUTPUT PARAM.).
243C              DEFAULT IS IWORK(4)=0.
244C
245C       THE FOLLOWING 3 PARAMETERS ARE IMPORTANT FOR
246C       DIFFERENTIAL-ALGEBRAIC SYSTEMS OF INDEX > 1.
247C       THE FUNCTION-SUBROUTINE SHOULD BE WRITTEN SUCH THAT
248C       THE INDEX 1,2,3 VARIABLES APPEAR IN THIS ORDER.
249C       IN ESTIMATING THE ERROR THE INDEX 2 VARIABLES ARE
250C       MULTIPLIED BY H, THE INDEX 3 VARIABLES BY H**2.
251C
252C    IWORK(5)  DIMENSION OF THE INDEX 1 VARIABLES (MUST BE > 0). FOR
253C              ODE'S THIS EQUALS THE DIMENSION OF THE SYSTEM.
254C              DEFAULT IWORK(5)=N.
255C
256C    IWORK(6)  DIMENSION OF THE INDEX 2 VARIABLES. DEFAULT IWORK(6)=0.
257C
258C    IWORK(7)  DIMENSION OF THE INDEX 3 VARIABLES. DEFAULT IWORK(7)=0.
259C
260C    IWORK(8)  SWITCH FOR STEP SIZE STRATEGY
261C              IF IWORK(8).EQ.1  MOD. PREDICTIVE CONTROLLER (GUSTAFSSON)
262C              IF IWORK(8).EQ.2  CLASSICAL STEP SIZE CONTROL
263C              THE DEFAULT VALUE (FOR IWORK(8)=0) IS IWORK(8)=1.
264C              THE CHOICE IWORK(8).EQ.1 SEEMS TO PRODUCE SAFER RESULTS;
265C              FOR SIMPLE PROBLEMS, THE CHOICE IWORK(8).EQ.2 PRODUCES
266C              OFTEN SLIGHTLY FASTER RUNS
267C
268C       IF THE DIFFERENTIAL SYSTEM HAS THE SPECIAL STRUCTURE THAT
269C            Y(I)' = Y(I+M2)   FOR  I=1,...,M1,
270C       WITH M1 A MULTIPLE OF M2, A SUBSTANTIAL GAIN IN COMPUTERTIME
271C       CAN BE ACHIEVED BY SETTING THE PARAMETERS IWORK(9) AND IWORK(10).
272C       E.G., FOR SECOND ORDER SYSTEMS P'=V, V'=G(P,V), WHERE P AND V ARE
273C       VECTORS OF DIMENSION N/2, ONE HAS TO PUT M1=M2=N/2.
274C       FOR M1>0 SOME OF THE INPUT PARAMETERS HAVE DIFFERENT MEANINGS:
275C       - JAC: ONLY THE ELEMENTS OF THE NON-TRIVIAL PART OF THE
276C              JACOBIAN HAVE TO BE STORED
277C              IF (MLJAC.EQ.N-M1) THE JACOBIAN IS SUPPOSED TO BE FULL
278C                 DFY(I,J) = PARTIAL F(I+M1) / PARTIAL Y(J)
279C                FOR I=1,N-M1 AND J=1,N.
280C              ELSE, THE JACOBIAN IS BANDED ( M1 = M2 * MM )
281C                 DFY(I-J+MUJAC+1,J+K*M2) = PARTIAL F(I+M1) / PARTIAL Y(J+K*M2)
282C                FOR I=1,MLJAC+MUJAC+1 AND J=1,M2 AND K=0,MM.
283C       - MLJAC: MLJAC=N-M1: IF THE NON-TRIVIAL PART OF THE JACOBIAN IS FULL
284C                0<=MLJAC<N-M1: IF THE (MM+1) SUBMATRICES (FOR K=0,MM)
285C                     PARTIAL F(I+M1) / PARTIAL Y(J+K*M2),  I,J=1,M2
286C                    ARE BANDED, MLJAC IS THE MAXIMAL LOWER BANDWIDTH
287C                    OF THESE MM+1 SUBMATRICES
288C       - MUJAC: MAXIMAL UPPER BANDWIDTH OF THESE MM+1 SUBMATRICES
289C                NEED NOT BE DEFINED IF MLJAC=N-M1
290C       - MAS: IF IMAS=0 THIS MATRIX IS ASSUMED TO BE THE IDENTITY AND
291C              NEED NOT BE DEFINED. SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
292C              IT IS ASSUMED THAT ONLY THE ELEMENTS OF RIGHT LOWER BLOCK OF
293C              DIMENSION N-M1 DIFFER FROM THAT OF THE IDENTITY MATRIX.
294C              IF (MLMAS.EQ.N-M1) THIS SUBMATRIX IS SUPPOSED TO BE FULL
295C                 AM(I,J) = M(I+M1,J+M1)     FOR I=1,N-M1 AND J=1,N-M1.
296C              ELSE, THE MASS MATRIX IS BANDED
297C                 AM(I-J+MUMAS+1,J) = M(I+M1,J+M1)
298C       - MLMAS: MLMAS=N-M1: IF THE NON-TRIVIAL PART OF M IS FULL
299C                0<=MLMAS<N-M1: LOWER BANDWIDTH OF THE MASS MATRIX
300C       - MUMAS: UPPER BANDWIDTH OF THE MASS MATRIX
301C                NEED NOT BE DEFINED IF MLMAS=N-M1
302C
303C    IWORK(9)  THE VALUE OF M1.  DEFAULT M1=0.
304C
305C    IWORK(10) THE VALUE OF M2.  DEFAULT M2=M1.
306C
307C    IWORK(11) NSMIN, MINIMAL NUMBER OF STAGES NS (ORDER 2*NS-1)
308C              POSSIBLE VALUES ARE 1,3,5,7. DEFAULT NS=3.
309C
310C    IWORK(12) NSMAX, MAXIMAL NUMBER OF STAGES NS.
311C              POSSIBLE VALUES ARE 1,3,5,7. DEFAULT NS=7.
312C
313C    IWORK(13) VALUE OF NS FOR THE FIRST STEP (DEFAULT VALUE: NSMIN)
314C
315C ----------
316C
317C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
318C
319C    WORK(2)   THE SAFETY FACTOR IN STEP SIZE PREDICTION,
320C              DEFAULT 0.9D0.
321C
322C    WORK(3)   DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
323C              INCREASE WORK(3), TO 0.1 SAY, WHEN JACOBIAN EVALUATIONS
324C              ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER
325C              (0.001D0, SAY). NEGATIV WORK(3) FORCES THE CODE TO
326C              COMPUTE THE JACOBIAN AFTER EVERY ACCEPTED STEP.
327C              DEFAULT 0.001D0.
328C
329C    WORK(5) AND WORK(6) : IF WORK(5) < HNEW/HOLD < WORK(6), THEN THE
330C              STEP SIZE IS NOT CHANGED. THIS SAVES, TOGETHER WITH A
331C              LARGE WORK(3), LU-DECOMPOSITIONS AND COMPUTING TIME FOR
332C              LARGE SYSTEMS. FOR SMALL SYSTEMS ONE MAY HAVE
333C              WORK(5)=1.D0, WORK(6)=1.2D0, FOR LARGE FULL SYSTEMS
334C              WORK(5)=0.99D0, WORK(6)=2.D0 MIGHT BE GOOD.
335C              DEFAULTS WORK(5)=1.D0, WORK(6)=1.2D0 .
336C
337C    WORK(7)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
338C
339C    WORK(8), WORK(9)   PARAMETERS FOR STEP SIZE SELECTION
340C              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION
341C                 WORK(8) <= HNEW/HOLD <= WORK(9)
342C              DEFAULT VALUES: WORK(8)=0.2D0, WORK(9)=8.D0
343C
344C    WORK(10)  ORDER IS INCREASED IF THE CONTRACTIVITY FACTOR IS
345C              SMALL THAN WORK(10), DEFAULT VALUE IS 0.002
346C
347C    WORK(11)  ORDER IS DECREASED IF THE CONTRACTIVITY FACTOR IS
348C              LARGER THAN WORK(11), DEFAULT VALUE IS 0.8
349C
350C    WORK(12), WORK(13)  ORDER IS DECREASED ONLY IF THE STEPSIZE
351C              RATIO SATISFIES  WORK(13).LE.HNEW/H.LE.WORK(12),
352C              DEFAULT VALUES ARE 1.2 AND 0.8
353C
354C-----------------------------------------------------------------------
355C
356C     OUTPUT PARAMETERS
357C     -----------------
358C     X           X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED
359C                 (AFTER SUCCESSFUL RETURN X=XEND).
360C
361C     Y(N)        NUMERICAL SOLUTION AT X
362C
363C     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
364C
365C     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
366C                   IDID= 1  COMPUTATION SUCCESSFUL,
367C                   IDID= 2  COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT)
368C                   IDID=-1  INPUT IS NOT CONSISTENT,
369C                   IDID=-2  LARGER NMAX IS NEEDED,
370C                   IDID=-3  STEP SIZE BECOMES TOO SMALL,
371C                   IDID=-4  MATRIX IS REPEATEDLY SINGULAR.
372C
373C   IWORK(14)  NFCN    NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL
374C                      EVALUATION OF THE JACOBIAN ARE NOT COUNTED)
375C   IWORK(15)  NJAC    NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
376C                      OR NUMERICALLY)
377C   IWORK(16)  NSTEP   NUMBER OF COMPUTED STEPS
378C   IWORK(17)  NACCPT  NUMBER OF ACCEPTED STEPS
379C   IWORK(18)  NREJCT  NUMBER OF REJECTED STEPS (DUE TO ERROR TEST),
380C                      (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED)
381C   IWORK(19)  NDEC    NUMBER OF LU-DECOMPOSITIONS OF THE MATRICES
382C   IWORK(20)  NSOL    NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS, OF ALL
383C                      SYSTEMS THAT HAVE TO BE SOLVED FOR ONE SIMPLIFIED
384C                      NEWTON ITERATION; THE NSTEP (REAL) FORWARD-BACKWARD
385C                      SUBSTITUTIONS, NEEDED FOR STEP SIZE SELECTION,
386C                      ARE NOT COUNTED.
387C-----------------------------------------------------------------------
388C *** *** *** *** *** *** *** *** *** *** *** *** ***
389C          DECLARATIONS
390C *** *** *** *** *** *** *** *** *** *** *** *** ***
391      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
392      DIMENSION Y(N),ATOL(*),RTOL(*),WORK(LWORK),IWORK(LIWORK)
393      DIMENSION RPAR(*),IPAR(*)
394      LOGICAL IMPLCT,JBAND,ARRET,STARTN,PRED
395      EXTERNAL FCN,JAC,MAS,SOLOUT
396C *** *** *** *** *** *** ***
397C        SETTING THE PARAMETERS
398C *** *** *** *** *** *** ***
399       NFCN=0
400       NJAC=0
401       NSTEP=0
402       NACCPT=0
403       NREJCT=0
404       NDEC=0
405       NSOL=0
406       ARRET=.FALSE.
407C -------- NUMBER MAXIMAL AND MINIMAL OF STAGES  NS
408      IF (IWORK(11).EQ.0) THEN
409         NSMIN=3
410      ELSE
411         NSMIN=MAX(1,IWORK(11))
412         IF (IWORK(11).GE.2) NSMIN=MAX(3,IWORK(11))
413         IF (IWORK(11).GE.4) NSMIN=MAX(5,IWORK(11))
414         IF (IWORK(11).GE.6) NSMIN=7
415      END IF
416      IF (IWORK(12).EQ.0) THEN
417         NSMAX=7
418      ELSE
419         NSMAX=MIN(7,IWORK(12))
420         IF (IWORK(12).LE.6) NSMAX=MIN(5,IWORK(12))
421         IF (IWORK(12).LE.4) NSMAX=MIN(3,IWORK(12))
422         IF (IWORK(12).LE.2) NSMAX=1
423      END IF
424      NS=NSMAX
425      IF (IWORK(13).EQ.0) THEN
426         NSUS=NSMIN
427      ELSE
428         NSUS=IWORK(13)
429         IF (NSUS.LE.0.OR.NS.GE.8.OR.NS.EQ.2.OR.NS.EQ.4.OR.NS.EQ.6) THEN
430            WRITE(6,*)' WRONG INPUT IWORK(13)=',IWORK(13)
431            ARRET=.TRUE.
432         END IF
433      END IF
434C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
435      IF (IWORK(2).EQ.0) THEN
436         NMAX=100000
437      ELSE
438         NMAX=IWORK(2)
439         IF (NMAX.LE.0) THEN
440            WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2)
441            ARRET=.TRUE.
442         END IF
443      END IF
444C -------- NIT    MAXIMAL NUMBER OF NEWTON ITERATIONS
445      IF (IWORK(3).EQ.0) THEN
446         NIT=7
447      ELSE
448         NIT=IWORK(3)
449         IF (NIT.LE.0.OR.NIT.GT.50) THEN
450            WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
451            ARRET=.TRUE.
452         END IF
453      END IF
454C -------- STARTN  SWITCH FOR STARTING VALUES OF NEWTON ITERATIONS
455      IF(IWORK(4).EQ.0)THEN
456         STARTN=.FALSE.
457      ELSE
458         STARTN=.TRUE.
459      END IF
460C -------- PARAMETER FOR DIFFERENTIAL-ALGEBRAIC COMPONENTS
461      NIND1=IWORK(5)
462      NIND2=IWORK(6)
463      NIND3=IWORK(7)
464      IF (NIND1.EQ.0) NIND1=N
465      IF (NIND1+NIND2+NIND3.NE.N) THEN
466       WRITE(6,*)' CURIOUS INPUT FOR IWORK(5,6,7)=',NIND1,NIND2,NIND3
467       ARRET=.TRUE.
468      END IF
469C -------- PRED   STEP SIZE CONTROL
470      IF(IWORK(8).LE.1)THEN
471         PRED=.TRUE.
472      ELSE
473         PRED=.FALSE.
474      END IF
475C -------- PARAMETER FOR SECOND ORDER EQUATIONS
476      M1=IWORK(9)
477      M2=IWORK(10)
478      NM1=N-M1
479      IF (M1.EQ.0) M2=N
480      IF (M2.EQ.0) M2=M1
481      IF (M1.LT.0.OR.M2.LT.0.OR.M1+M2.GT.N) THEN
482       WRITE(6,*)' CURIOUS INPUT FOR IWORK(9,10)=',M1,M2
483       ARRET=.TRUE.
484      END IF
485C -------- UROUND   SMALLEST NUMBER SATISFYING 1.0D0+UROUND>1.0D0
486      IF (WORK(1).EQ.0.0D0) THEN
487         UROUND=1.0D-16
488      ELSE
489         UROUND=WORK(1)
490         IF (UROUND.LE.1.0D-19.OR.UROUND.GE.1.0D0) THEN
491            WRITE(6,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',WORK(1)
492            ARRET=.TRUE.
493         END IF
494      END IF
495C --------- CHECK IF TOLERANCES ARE O.K.
496      IF (ITOL.EQ.0) THEN
497          IF (ATOL(1).LE.0.D0.OR.RTOL(1).LE.10.D0*UROUND) THEN
498              WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
499              ARRET=.TRUE.
500          END IF
501      ELSE
502          DO I=1,N
503          IF (ATOL(I).LE.0.D0.OR.RTOL(I).LE.10.D0*UROUND) THEN
504              WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
505              ARRET=.TRUE.
506          END IF
507          END DO
508      END IF
509C --------- SAFE     SAFETY FACTOR IN STEP SIZE PREDICTION
510      IF (WORK(2).EQ.0.0D0) THEN
511         SAFE=0.9D0
512      ELSE
513         SAFE=WORK(2)
514         IF (SAFE.LE.0.001D0.OR.SAFE.GE.1.0D0) THEN
515            WRITE(6,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2)
516            ARRET=.TRUE.
517         END IF
518      END IF
519C ------ THET     DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
520      IF (WORK(3).EQ.0.D0) THEN
521         THET=0.001D0
522      ELSE
523         THET=WORK(3)
524         IF (THET.GE.1.0D0) THEN
525            WRITE(6,*)' CURIOUS INPUT FOR WORK(3)=',WORK(3)
526            ARRET=.TRUE.
527         END IF
528      END IF
529C --- QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST.
530      IF (WORK(5).EQ.0.D0) THEN
531         QUOT1=1.D0
532      ELSE
533         QUOT1=WORK(5)
534      END IF
535      IF (WORK(6).EQ.0.D0) THEN
536         QUOT2=1.2D0
537      ELSE
538         QUOT2=WORK(6)
539      END IF
540      IF (QUOT1.GT.1.0D0.OR.QUOT2.LT.1.0D0) THEN
541         WRITE(6,*)' CURIOUS INPUT FOR WORK(5,6)=',QUOT1,QUOT2
542         ARRET=.TRUE.
543      END IF
544C -------- MAXIMAL STEP SIZE
545      IF (WORK(7).EQ.0.D0) THEN
546         HMAX=XEND-X
547      ELSE
548         HMAX=WORK(7)
549      END IF
550C -------  FACL,FACR     PARAMETERS FOR STEP SIZE SELECTION
551      IF(WORK(8).EQ.0.D0)THEN
552         FACL=5.D0
553      ELSE
554         FACL=1.D0/WORK(8)
555      END IF
556      IF(WORK(9).EQ.0.D0)THEN
557         FACR=1.D0/8.0D0
558      ELSE
559         FACR=1.D0/WORK(9)
560      END IF
561      IF (FACL.LT.1.0D0.OR.FACR.GT.1.0D0) THEN
562            WRITE(6,*)' CURIOUS INPUT WORK(8,9)=',WORK(8),WORK(9)
563            ARRET=.TRUE.
564         END IF
565C -------- PARAMETERS FOR ORDER SELECTION STRATEGY
566      IF (WORK(10).EQ.0.D0) THEN
567         VITU=0.002D0
568      ELSE
569         VITU=WORK(10)
570      END IF
571      IF (WORK(11).EQ.0.D0) THEN
572         VITD=0.8D0
573      ELSE
574         VITD=WORK(11)
575      END IF
576      IF (WORK(12).EQ.0.D0) THEN
577         HHOU=1.2D0
578      ELSE
579         HHOU=WORK(12)
580      END IF
581      IF (WORK(13).EQ.0.D0) THEN
582         HHOD=0.8D0
583      ELSE
584         HHOD=WORK(13)
585      END IF
586C *** *** *** *** *** *** *** *** *** *** *** *** ***
587C         COMPUTATION OF ARRAY ENTRIES
588C *** *** *** *** *** *** *** *** *** *** *** *** ***
589C ---- IMPLICIT, BANDED OR NOT ?
590      IMPLCT=IMAS.NE.0
591      JBAND=MLJAC.LT.NM1
592C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
593C -- JACOBIAN  AND  MATRICES E1, E2
594      IF (JBAND) THEN
595         LDJAC=MLJAC+MUJAC+1
596         LDE1=MLJAC+LDJAC
597      ELSE
598         MLJAC=NM1
599         MUJAC=NM1
600         LDJAC=NM1
601         LDE1=NM1
602      END IF
603C -- MASS MATRIX
604      IF (IMPLCT) THEN
605          IF (MLMAS.NE.NM1) THEN
606              LDMAS=MLMAS+MUMAS+1
607              IF (JBAND) THEN
608                 IJOB=4
609              ELSE
610                 IJOB=3
611              END IF
612          ELSE
613              LDMAS=NM1
614              IJOB=5
615          END IF
616C ------ BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF "JAC"
617          IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN
618             WRITE (6,*) 'BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF
619     & "JAC"'
620            ARRET=.TRUE.
621          END IF
622      ELSE
623          LDMAS=0
624          IF (JBAND) THEN
625             IJOB=2
626          ELSE
627             IJOB=1
628             IF (N.GT.2.AND.IWORK(1).NE.0) IJOB=7
629          END IF
630      END IF
631      LDMAS2=MAX(1,LDMAS)
632C ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN
633      IF ((IMPLCT.OR.JBAND).AND.IJOB.EQ.7) THEN
634         WRITE(6,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS WITH
635     &FULL JACOBIAN'
636         ARRET=.TRUE.
637      END IF
638C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
639      NNS=NS*N
640      NM1NS=NS*NM1
641      NMEE=(NS-1)*NM1
642      IEZZ=21
643      IEY0=IEZZ+NNS
644      IESCAL=IEY0+N
645      IEFF=IESCAL+N
646      IECON=IEFF+NNS
647      IEJAC=IECON+NNS+N
648      IEMAS=IEJAC+N*LDJAC
649      IEE1=IEMAS+NM1*LDMAS
650      IEE=IEE1+NM1*LDE1
651C ------ TOTAL STORAGE REQUIREMENT -----------
652      ISTORE=IEE+NMEE*LDE1-1
653      IF(ISTORE.GT.LWORK)THEN
654         WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
655         ARRET=.TRUE.
656      END IF
657C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
658      IEIP1=21
659      IEIP2=IEIP1+NM1
660      IEIPH=IEIP2+NM1*(NS-1)/2
661C --------- TOTAL REQUIREMENT ---------------
662      ISTORE=IEIPH+NM1-1
663      IF (ISTORE.GT.LIWORK) THEN
664         WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
665         ARRET=.TRUE.
666      END IF
667C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
668      IF (ARRET) THEN
669         IDID=-1
670         RETURN
671      END IF
672C -------- CALL TO CORE INTEGRATOR ------------
673      CALL RADCOV(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,NSUS,
674     &   JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
675     &   NMAX,UROUND,SAFE,THET,QUOT1,QUOT2,NIT,IJOB,STARTN,
676     &   NIND1,NIND2,NIND3,PRED,FACL,FACR,M1,M2,NM1,NSMIN,NS,NNS,NM1NS,
677     &   NMEE,IMPLCT,JBAND,LDJAC,LDE1,LDMAS2,WORK(IEZZ),WORK(IEY0),
678     &   WORK(IESCAL),WORK(IEFF),WORK(IEJAC),WORK(IEE1),WORK(IEE),
679     &   WORK(IEMAS),WORK(IECON),IWORK(IEIP1),IWORK(IEIP2),
680     &   IWORK(IEIPH),VITU,VITD,HHOU,HHOD,NFCN,NJAC,NSTEP,NACCPT,
681     &   NREJCT,NDEC,NSOL,RPAR,IPAR)
682      IWORK(13)=NSUS
683      IWORK(14)=NFCN
684      IWORK(15)=NJAC
685      IWORK(16)=NSTEP
686      IWORK(17)=NACCPT
687      IWORK(18)=NREJCT
688      IWORK(19)=NDEC
689      IWORK(20)=NSOL
690C ----------- RETURN -----------
691      RETURN
692      END
693C
694C     END OF SUBROUTINE RADAU
695C
696C ***********************************************************
697C
698      SUBROUTINE RADCOV(N,FCN,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,NS,
699     &   JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
700     &   NMAX,UROUND,SAFE,THET,QUOT1,QUOT2,NIT1,IJOB,STARTN,NIND1,
701     &   NIND2,NIND3,PRED,FACL,FACR,M1,M2,NM1,NSMIN,NSMAX,NNMS,NM1NS,
702     &   NMEE,IMPLCT,BANDED,LDJAC,LDE1,LDMAS,ZZ,Y0,SCAL,FF,FJAC,
703     &   E1,EE2,FMAS,CONT,IP1,IP2,IPHES,VITU,VITD,HHOU,HHOD,
704     &   NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,RPAR,IPAR)
705C ----------------------------------------------------------
706C     CORE INTEGRATOR FOR RADAU
707C     PARAMETERS SAME AS IN RADAU WITH WORKSPACE ADDED
708C ----------------------------------------------------------
709C         DECLARATIONS
710C ----------------------------------------------------------
711      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
712      PARAMETER (NSDIM=7)
713C --- THIS PARAMETER HAS TO BE CHANGED IF NUMBER OF STAGES IS >=7
714      DIMENSION Y(N),ZZ(NNMS),FF(NNMS),Y0(N),SCAL(N)
715      DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),CONT(NNMS+N)
716      DIMENSION ATOL(1),RTOL(1),RPAR(*),IPAR(*)
717      DIMENSION E1(LDE1,NM1),EE2(LDE1,NMEE)
718      DIMENSION ALPH(NSDIM),BETA(NSDIM),DD(NSDIM)
719      DIMENSION ALPHN(NSDIM),BETAN(NSDIM)
720      INTEGER IP1(NM1),IP2(NMEE/2),IPHES(NM1)
721      COMMON/WEIGHT/NN,NSCON,XSOL,HSOL,C(0:NSDIM)
722      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
723      COMMON /COE3/ T311,T312,T313,T321,T322,T323,T331,
724     *    TI311,TI312,TI313,TI321,TI322,TI323,TI331,TI332,TI333
725      COMMON /COE5/ T511,T512,T513,T514,T515,T521,T522,T523,T524,T525,
726     *    T531,T532,T533,T534,T535,T541,T542,T543,T544,T545,T551,
727     *    TI511,TI512,TI513,TI514,TI515,TI521,TI522,TI523,TI524,TI525,
728     *    TI531,TI532,TI533,TI534,TI535,TI541,TI542,TI543,TI544,TI545,
729     *    TI551,TI552,TI553,TI554,TI555
730      COMMON /COE7/ T711,T712,T713,T714,T715,T716,T717,
731     *              T721,T722,T723,T724,T725,T726,T727,
732     *              T731,T732,T733,T734,T735,T736,T737,
733     *              T741,T742,T743,T744,T745,T746,T747,
734     *              T751,T752,T753,T754,T755,T756,T757,
735     *              T761,T762,T763,T764,T765,T766,T767,T771,
736     *              TI711,TI712,TI713,TI714,TI715,TI716,TI717,
737     *              TI721,TI722,TI723,TI724,TI725,TI726,TI727,
738     *              TI731,TI732,TI733,TI734,TI735,TI736,TI737,
739     *              TI741,TI742,TI743,TI744,TI745,TI746,TI747,
740     *              TI751,TI752,TI753,TI754,TI755,TI756,TI757,
741     *              TI761,TI762,TI763,TI764,TI765,TI766,TI767,
742     *              TI771,TI772,TI773,TI774,TI775,TI776,TI777
743      LOGICAL REJECT,FIRST,IMPLCT,BANDED,CALJAC,STARTN,CALHES
744      LOGICAL INDEX1,INDEX2,INDEX3,LAST,PRED,CHANGE,UNEXP,UNEXN,VARIAB
745      EXTERNAL FCN,JAC
746C *** *** *** *** *** *** ***
747C  INITIALISATIONS
748C *** *** *** *** *** *** ***
749C -------- CHECK THE INDEX OF THE PROBLEM -----
750      INDEX1=NIND1.NE.0
751      INDEX2=NIND2.NE.0
752      INDEX3=NIND3.NE.0
753C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
754      IF (IMPLCT) CALL MAS(NM1,FMAS,LDMAS,RPAR,IPAR)
755      VARIAB=NSMIN.LT.NSMAX
756C ---------- CONSTANTS ---------
757      EXPO=1.D0/(NS+1.D0)
758      SQ6=DSQRT(6.D0)
759      C31=(4.D0-SQ6)/10.D0
760      C32=(4.D0+SQ6)/10.D0
761      C31M1=C31-1.D0
762      C32M1=C32-1.D0
763      C31MC2=C31-C32
764      DD1=-(13.D0+7.D0*SQ6)/3.D0
765      DD2=(-13.D0+7.D0*SQ6)/3.D0
766      DD3=-1.D0/3.D0
767      N2=2*N
768      N3=3*N
769      N4=4*N
770      N5=5*N
771      N6=6*N
772      UNEXP=.FALSE.
773      UNEXN=.FALSE.
774      CHANGE=.FALSE.
775      IKEEP=0
776      ICHAN=0
777      THETA=0.0D0
778      THETAT=0.0D0
779      NN=N
780      NNS=N*NS
781      NSCON=NS
782      LRC=NNS+N
783      CALL COERTV(NSMAX)
784      CALL COERCV(NS,C,DD,U1,ALPH,BETA)
785      IF (M1.GT.0) IJOB=IJOB+10
786      POSNEG=SIGN(1.D0,XEND-X)
787      HMAXN=MIN(ABS(HMAX),ABS(XEND-X))
788      IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6
789      H=MIN(ABS(H),HMAXN)
790      H=SIGN(H,POSNEG)
791      HOLD=H
792      REJECT=.FALSE.
793      FIRST=.TRUE.
794      LAST=.FALSE.
795      IF ((X+H*1.0001D0-XEND)*POSNEG.GE.0.D0) THEN
796         H=XEND-X
797         LAST=.TRUE.
798      END IF
799      HOPT=H
800      FACCON=1.D0
801      NSING=0
802      XOLD=X
803      IF (IOUT.NE.0) THEN
804          IRTRN=1
805          NRSOL=1
806          XOSOL=XOLD
807          XSOL=X
808          DO I=1,N
809             CONT(I)=Y(I)
810          END DO
811          NSOLU=N
812          HSOL=HOLD
813          CALL SOLOUT(NRSOL,XOSOL,XSOL,Y,CONT,LRC,NSOLU,
814     &                RPAR,IPAR,IRTRN)
815          IF (IRTRN.LT.0) GOTO 179
816      END IF
817      MLE=MLJAC
818      MUE=MUJAC
819      MBJAC=MLJAC+MUJAC+1
820      MBB=MLMAS+MUMAS+1
821      MDIAG=MLE+MUE+1
822      MDIFF=MLE+MUE-MUMAS
823      MBDIAG=MUMAS+1
824      EXPMNS=(NS+1.0D0)/(2.0D0*NS)
825      QUOTT=ATOL(1)/RTOL(1)
826      IF (ITOL.EQ.0) THEN
827          RTOL1=0.1D0*RTOL(1)**EXPMNS
828          ATOL1=RTOL1*QUOTT
829          DO I=1,N
830             SCAL(I)=ATOL1+RTOL1*ABS(Y(I))
831          END DO
832      ELSE
833          DO I=1,N
834             QUOTT=ATOL(I)/RTOL(I)
835             RTOL1=0.1D0*RTOL(I)**EXPMNS
836             ATOL1=RTOL1*QUOTT
837             SCAL(I)=ATOL1+RTOL1*ABS(Y(I))
838          END DO
839      END IF
840      HHFAC=H
841      CALL FCN(N,X,Y,Y0,RPAR,IPAR)
842      NFCN=NFCN+1
843C --- BASIC INTEGRATION STEP
844  10  CONTINUE
845C *** *** *** *** *** *** ***
846C  COMPUTATION OF THE JACOBIAN
847C *** *** *** *** *** *** ***
848      NJAC=NJAC+1
849      IF (IJAC.EQ.0) THEN
850C --- COMPUTE JACOBIAN MATRIX NUMERICALLY
851         IF (BANDED) THEN
852C --- JACOBIAN IS BANDED
853            MUJACP=MUJAC+1
854            MD=MIN(MBJAC,M2)
855            DO MM=1,M1/M2+1
856               DO K=1,MD
857                  J=K+(MM-1)*M2
858 12               FF(J)=Y(J)
859                  FF(J+N)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J))))
860                  Y(J)=Y(J)+FF(J+N)
861                  J=J+MD
862                  IF (J.LE.MM*M2) GOTO 12
863                  CALL FCN(N,X,Y,CONT,RPAR,IPAR)
864                  J=K+(MM-1)*M2
865                  J1=K
866                  LBEG=MAX(1,J1-MUJAC)+M1
867 14               LEND=MIN(M2,J1+MLJAC)+M1
868                  Y(J)=FF(J)
869                  MUJACJ=MUJACP-J1-M1
870                  DO L=LBEG,LEND
871                     FJAC(L+MUJACJ,J)=(CONT(L)-Y0(L))/FF(J+N)
872                  END DO
873                  J=J+MD
874                  J1=J1+MD
875                  LBEG=LEND+1
876                  IF (J.LE.MM*M2) GOTO 14
877               END DO
878            END DO
879         ELSE
880C --- JACOBIAN IS FULL
881            DO I=1,N
882               YSAFE=Y(I)
883               DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE)))
884               Y(I)=YSAFE+DELT
885               CALL FCN(N,X,Y,CONT,RPAR,IPAR)
886               DO J=M1+1,N
887                 FJAC(J-M1,I)=(CONT(J)-Y0(J))/DELT
888               END DO
889               Y(I)=YSAFE
890            END DO
891         END IF
892      ELSE
893C --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
894         CALL JAC(N,X,Y,FJAC,LDJAC,RPAR,IPAR)
895      END IF
896      CALJAC=.TRUE.
897      CALHES=.TRUE.
898  20  CONTINUE
899C --- CHANGE THE ORDER HERE IF NECESSARY
900      IF (VARIAB) THEN
901         NSNEW=NS
902         ICHAN=ICHAN+1
903         HQUOT=H/HOLD
904         THETAT=MIN(10.0D0,MAX(THETA,THETAT*0.5))
905         IF (NEWT.GT.1.AND.THETAT.LE.VITU.AND.
906     &      HQUOT.LT.HHOU.AND.HQUOT.GT.HHOD) NSNEW=MIN(NSMAX,NS+2)
907         IF (THETAT.GE.VITD.OR.UNEXP) NSNEW=MAX(NSMIN,NS-2)
908         IF (ICHAN.GE.1.AND.UNEXN) NSNEW=MAX(NSMIN,NS-2)
909         IF (ICHAN.LE.10) NSNEW=MIN(NS,NSNEW)
910         CHANGE=NS.NE.NSNEW
911         UNEXN=.FALSE.
912         UNEXP=.FALSE.
913         IF (CHANGE) THEN
914            NS=NSNEW
915            ICHAN=1
916            NNS=N*NS
917            NSCON=NS
918            LRC=NNS+N
919            CALL COERCV(NS,C,DD,U1,ALPH,BETA)
920            EXPO=1.D0/(NS+1.D0)
921            EXPMNS=(NS+1.0D0)/(2.0D0*NS)
922            RTOL1=0.1D0*RTOL(1)**EXPMNS
923            ATOL1=RTOL1*QUOTT
924            IF (ITOL.EQ.0) THEN
925               DO I=1,N
926                  SCAL(I)=ATOL1+RTOL1*ABS(Y(I))
927               END DO
928            ELSE
929               DO I=1,N
930                  QUOTT=ATOL(I)/RTOL(I)
931                  RTOL1=0.1D0*RTOL(I)**EXPMNS
932                  ATOL1=RTOL1*QUOTT
933                  SCAL(I)=ATOL1+RTOL1*ABS(Y(I))
934               END DO
935            END IF
936         END IF
937      END IF
938C --- COMPUTE THE MATRICES E1 AND E2 AND THEIR DECOMPOSITIONS
939      FAC1=U1/H
940      CALL DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
941     &            M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES)
942      IF (IER.NE.0) GOTO 78
943      DO K=1,(NS-1)/2
944         ALPHN(K)=ALPH(K)/H
945         BETAN(K)=BETA(K)/H
946         IAD=(K-1)*2*NM1+1
947         CALL DECOMC(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
948     &         M1,M2,NM1,ALPHN(K),BETAN(K),EE2(1,IAD),EE2(1,IAD+NM1),
949     &         LDE1,IP2((K-1)*NM1+1),IER,IJOB)
950         IF (IER.NE.0) GOTO 78
951      END DO
952      NDEC=NDEC+1
953  30  CONTINUE
954      IF (VARIAB.AND.IKEEP.EQ.1) THEN
955         ICHAN=ICHAN+1
956         IKEEP=0
957         IF (ICHAN.GE.10.AND.NS.LT.NSMAX) GOTO 20
958      END IF
959      NSTEP=NSTEP+1
960      IF (NSTEP.GT.NMAX) GOTO 178
961      IF (0.1D0*ABS(H).LE.ABS(X)*UROUND) GOTO 177
962          IF (INDEX2) THEN
963             DO I=NIND1+1,NIND1+NIND2
964                SCAL(I)=SCAL(I)/HHFAC
965             END DO
966          END IF
967          IF (INDEX3) THEN
968             DO I=NIND1+NIND2+1,NIND1+NIND2+NIND3
969                SCAL(I)=SCAL(I)/(HHFAC*HHFAC)
970             END DO
971          END IF
972      XPH=X+H
973C *** *** *** *** *** *** ***
974C *** *** *** *** *** *** ***
975      IF (NS.EQ.3) THEN
976C *** *** *** *** *** *** ***
977C *** *** *** *** *** *** ***
978      IF (FIRST.OR.STARTN.OR.CHANGE) THEN
979         DO I=1,NNS
980            ZZ(I)=0.D0
981            FF(I)=0.D0
982         END DO
983      ELSE
984         HQUOT=H/HOLD
985         C3Q=HQUOT
986         C1Q=C31*C3Q
987         C2Q=C32*C3Q
988         DO I=1,N
989            AK1=CONT(I+N)
990            AK2=CONT(I+N2)
991            AK3=CONT(I+N3)
992            Z1I=C1Q*(AK1+(C1Q-C32M1)*(AK2+(C1Q-C31M1)*AK3))
993            Z2I=C2Q*(AK1+(C2Q-C32M1)*(AK2+(C2Q-C31M1)*AK3))
994            Z3I=C3Q*(AK1+(C3Q-C32M1)*(AK2+(C3Q-C31M1)*AK3))
995            ZZ(I)=Z1I
996            ZZ(I+N)=Z2I
997            ZZ(I+N2)=Z3I
998            FF(I)=TI311*Z1I+TI312*Z2I+TI313*Z3I
999            FF(I+N)=TI321*Z1I+TI322*Z2I+TI323*Z3I
1000            FF(I+N2)=TI331*Z1I+TI332*Z2I+TI333*Z3I
1001         END DO
1002      END IF
1003C *** *** *** *** *** *** ***
1004C  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
1005C *** *** *** *** *** *** ***
1006      NEWT=0
1007      NIT=NIT1
1008      EXPMI=1.0D0/EXPMNS
1009      FNEWT=MAX(10*UROUND/RTOL1,MIN(0.03D0,RTOL1**(EXPMI-1.0D0)))
1010      FACCON=MAX(FACCON,UROUND)**0.8D0
1011      THETA=ABS(THET)
1012  40  CONTINUE
1013      IF (NEWT.GE.NIT) GOTO 78
1014C ---     COMPUTE THE RIGHT-HAND SIDE
1015      DO I=1,N
1016         CONT(I)=Y(I)+ZZ(I)
1017      END DO
1018      CALL FCN(N,X+C31*H,CONT,ZZ,RPAR,IPAR)
1019      DO I=1,N
1020         CONT(I)=Y(I)+ZZ(I+N)
1021      END DO
1022      CALL FCN(N,X+C32*H,CONT,ZZ(1+N),RPAR,IPAR)
1023      DO I=1,N
1024         CONT(I)=Y(I)+ZZ(I+N2)
1025      END DO
1026      CALL FCN(N,XPH,CONT,ZZ(1+N2),RPAR,IPAR)
1027      NFCN=NFCN+3
1028C ---     SOLVE THE LINEAR SYSTEMS
1029      DO I=1,N
1030         A1=ZZ(I)
1031         A2=ZZ(I+N)
1032         A3=ZZ(I+N2)
1033         ZZ(I)=TI311*A1+TI312*A2+TI313*A3
1034         ZZ(I+N)=TI321*A1+TI322*A2+TI323*A3
1035         ZZ(I+N2)=TI331*A1+TI332*A2+TI333*A3
1036      END DO
1037      CALL SLVRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1038     &      M1,M2,NM1,FAC1,ALPHN(1),BETAN(1),E1,EE2,EE2(1,1+NM1),LDE1,
1039     &          ZZ,ZZ(1+N),ZZ(1+N2),FF,FF(1+N),FF(1+N2),
1040     &          CONT,IP1,IP2,IPHES,IER,IJOB)
1041      NSOL=NSOL+1
1042      NEWT=NEWT+1
1043      DYNO=0.D0
1044      DO I=1,N
1045         DENOM=SCAL(I)
1046         DYNO=DYNO+(ZZ(I)/DENOM)**2+(ZZ(I+N)/DENOM)**2
1047     &          +(ZZ(I+N2)/DENOM)**2
1048      END DO
1049      DYNO=DSQRT(DYNO/NNS)
1050C ---     BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE
1051      IF (NEWT.GT.1.AND.NEWT.LT.NIT) THEN
1052         THQ=DYNO/DYNOLD
1053         IF (NEWT.EQ.2) THEN
1054            THETA=THQ
1055         ELSE
1056            THETA=SQRT(THQ*THQOLD)
1057         END IF
1058         THQOLD=THQ
1059         IF (THETA.LT.0.99D0) THEN
1060            FACCON=THETA/(1.0D0-THETA)
1061            DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)/FNEWT
1062            IF (DYTH.GE.1.0D0) THEN
1063               QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH))
1064               HHFAC=0.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT))
1065               H=HHFAC*H
1066               REJECT=.TRUE.
1067               LAST=.FALSE.
1068               IF (HHFAC.LE.0.5D0) UNEXN=.TRUE.
1069               IF (CALJAC) GOTO 20
1070               GOTO 10
1071            END IF
1072         ELSE
1073            GOTO 78
1074         END IF
1075      END IF
1076      DYNOLD=MAX(DYNO,UROUND)
1077      DO I=1,N
1078         IN=I+N
1079         IN2=I+N2
1080         F1I=FF(I)+ZZ(I)
1081         F2I=FF(IN)+ZZ(IN)
1082         F3I=FF(IN2)+ZZ(IN2)
1083         FF(I)=F1I
1084         FF(IN)=F2I
1085         FF(IN2)=F3I
1086         ZZ(I)=T311*F1I+T312*F2I+T313*F3I
1087         ZZ(IN)=T321*F1I+T322*F2I+T323*F3I
1088         ZZ(IN2)=T331*F1I+    F2I
1089      END DO
1090      IF (FACCON*DYNO.GT.FNEWT) GOTO 40
1091C --- ERROR ESTIMATION
1092      CALL ESTRAD (N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1093     &          H,DD1,DD2,DD3,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,
1094     &          E1,LDE1,ZZ,ZZ(1+N),ZZ(1+N2),CONT,FF,FF(1+N),IP1,
1095     &          IPHES,SCAL,ERR,FIRST,REJECT,FAC1,RPAR,IPAR)
1096C       --- COMPUTE FINITE DIFFERENCES FOR DENSE OUTPUT
1097      IF (ERR.LT.1.D0) THEN
1098         DO I=1,N
1099            Y(I)=Y(I)+ZZ(I+N2)
1100            Z2I=ZZ(I+N)
1101            Z1I=ZZ(I)
1102            CONT(I+N)=(Z2I-ZZ(I+N2))/C32M1
1103            AK=(Z1I-Z2I)/C31MC2
1104            ACONT3=Z1I/C31
1105            ACONT3=(AK-ACONT3)/C32
1106            CONT(I+N2)=(AK-CONT(I+N))/C31M1
1107            CONT(I+N3)=CONT(I+N2)-ACONT3
1108         END DO
1109      END IF
1110C *** *** *** *** *** *** ***
1111C *** *** *** *** *** *** ***
1112      ELSE
1113      IF (NS.EQ.5) THEN
1114C *** *** *** *** *** *** ***
1115C *** *** *** *** *** *** ***
1116      IF (FIRST.OR.STARTN.OR.CHANGE) THEN
1117         DO I=1,NNS
1118            ZZ(I)=0.D0
1119            FF(I)=0.D0
1120         END DO
1121      ELSE
1122         HQUOT=H/HOLD
1123         DO K=1,NS
1124            CCQ=C(K)*HQUOT
1125            DO I=1,N
1126               VAL=CONT(I+NS*N)
1127               DO L=NS-1,1,-1
1128                  VAL=CONT(I+L*N)+(CCQ-C(NS-L)+1.D0)*VAL
1129               END DO
1130               ZZ(I+(K-1)*N)=CCQ*VAL
1131            END DO
1132         END DO
1133         DO I=1,N
1134            Z1I=ZZ(I)
1135            Z2I=ZZ(I+N)
1136            Z3I=ZZ(I+N2)
1137            Z4I=ZZ(I+N3)
1138            Z5I=ZZ(I+N4)
1139            FF(I)=TI511*Z1I+TI512*Z2I+TI513*Z3I+TI514*Z4I+TI515*Z5I
1140            FF(I+N)=TI521*Z1I+TI522*Z2I+TI523*Z3I+TI524*Z4I+TI525*Z5I
1141            FF(I+N2)=TI531*Z1I+TI532*Z2I+TI533*Z3I+TI534*Z4I+TI535*Z5I
1142            FF(I+N3)=TI541*Z1I+TI542*Z2I+TI543*Z3I+TI544*Z4I+TI545*Z5I
1143            FF(I+N4)=TI551*Z1I+TI552*Z2I+TI553*Z3I+TI554*Z4I+TI555*Z5I
1144         END DO
1145      END IF
1146C *** *** *** *** *** *** ***
1147C  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
1148C *** *** *** *** *** *** ***
1149      NEWT=0
1150      NIT=NIT1+5
1151      EXPMI=1.0D0/EXPMNS
1152      FNEWT=MAX(10*UROUND/RTOL1,MIN(0.03D0,RTOL1**(EXPMI-1.0D0)))
1153      FACCON=MAX(FACCON,UROUND)**0.8D0
1154      THETA=ABS(THET)
1155 140  CONTINUE
1156      IF (NEWT.GE.NIT) GOTO 78
1157C ---     COMPUTE THE RIGHT-HAND SIDE
1158      DO K=0,NS-1
1159         IADD=K*N
1160         DO I=1,N
1161            CONT(I)=Y(I)+ZZ(IADD+I)
1162         END DO
1163         CALL FCN(N,X+C(K+1)*H,CONT,ZZ(IADD+1),RPAR,IPAR)
1164      END DO
1165      NFCN=NFCN+NS
1166C ---     SOLVE THE LINEAR SYSTEMS
1167      DO I=1,N
1168         Z1I=ZZ(I)
1169         Z2I=ZZ(I+N)
1170         Z3I=ZZ(I+N2)
1171         Z4I=ZZ(I+N3)
1172         Z5I=ZZ(I+N4)
1173         ZZ(I)=TI511*Z1I+TI512*Z2I+TI513*Z3I+TI514*Z4I+TI515*Z5I
1174         ZZ(I+N)=TI521*Z1I+TI522*Z2I+TI523*Z3I+TI524*Z4I+TI525*Z5I
1175         ZZ(I+N2)=TI531*Z1I+TI532*Z2I+TI533*Z3I+TI534*Z4I+TI535*Z5I
1176         ZZ(I+N3)=TI541*Z1I+TI542*Z2I+TI543*Z3I+TI544*Z4I+TI545*Z5I
1177         ZZ(I+N4)=TI551*Z1I+TI552*Z2I+TI553*Z3I+TI554*Z4I+TI555*Z5I
1178      END DO
1179      CALL SLVRAR(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1180     &          M1,M2,NM1,FAC1,E1,LDE1,ZZ,FF,IP1,IPHES,IER,IJOB)
1181      DO K=1,2
1182         IAD=(K-1)*2*NM1+1
1183         CALL SLVRAI(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1184     &     M1,M2,NM1,ALPHN(K),BETAN(K),EE2(1,IAD),EE2(1,IAD+NM1),LDE1,
1185     &     ZZ(1+(2*K-1)*N),ZZ(1+2*K*N),FF(1+(2*K-1)*N),
1186     &     FF(1+2*K*N),CONT,IP2((K-1)*NM1+1),IPHES,IER,IJOB)
1187      END DO
1188      NSOL=NSOL+1
1189      NEWT=NEWT+1
1190      DYNO=0.D0
1191      DO I=1,N
1192         DENOM=SCAL(I)
1193         DO K=0,NS-1
1194            DYNO=DYNO+(ZZ(I+K*N)/DENOM)**2
1195         END DO
1196      END DO
1197      DYNO=DSQRT(DYNO/NNS)
1198C ---     BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE
1199      IF (NEWT.GT.1.AND.NEWT.LT.NIT) THEN
1200         THQ=DYNO/DYNOLD
1201         IF (NEWT.EQ.2) THEN
1202            THETA=THQ
1203         ELSE
1204            THETA=SQRT(THQ*THQOLD)
1205         END IF
1206         THQOLD=THQ
1207         IF (THETA.LT.0.99D0) THEN
1208            FACCON=THETA/(1.0D0-THETA)
1209            DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)/FNEWT
1210            IF (DYTH.GE.1.0D0) THEN
1211               QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH))
1212               HHFAC=0.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT))
1213               H=HHFAC*H
1214               REJECT=.TRUE.
1215               LAST=.FALSE.
1216               IF (HHFAC.LE.0.5D0) UNEXN=.TRUE.
1217               IF (CALJAC) GOTO 20
1218               GOTO 10
1219            END IF
1220         ELSE
1221            GOTO 78
1222         END IF
1223      END IF
1224      DYNOLD=MAX(DYNO,UROUND)
1225      DO I=1,N
1226         Z1I=FF(I)+ZZ(I)
1227         Z2I=FF(I+N)+ZZ(I+N)
1228         Z3I=FF(I+N2)+ZZ(I+N2)
1229         Z4I=FF(I+N3)+ZZ(I+N3)
1230         Z5I=FF(I+N4)+ZZ(I+N4)
1231         FF(I)=Z1I
1232         FF(I+N)=Z2I
1233         FF(I+N2)=Z3I
1234         FF(I+N3)=Z4I
1235         FF(I+N4)=Z5I
1236         ZZ(I)=T511*Z1I+T512*Z2I+T513*Z3I+T514*Z4I+T515*Z5I
1237         ZZ(I+N)=T521*Z1I+T522*Z2I+T523*Z3I+T524*Z4I+T525*Z5I
1238         ZZ(I+N2)=T531*Z1I+T532*Z2I+T533*Z3I+T534*Z4I+T535*Z5I
1239         ZZ(I+N3)=T541*Z1I+T542*Z2I+T543*Z3I+T544*Z4I+T545*Z5I
1240         ZZ(I+N4)=T551*Z1I+     Z2I         +     Z4I
1241      END DO
1242      IF (FACCON*DYNO.GT.FNEWT) GOTO 140
1243C --- ERROR ESTIMATION
1244      CALL ESTRAV (N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1245     &          H,DD,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,NS,NNS,E1,LDE1,
1246     &          ZZ,CONT,FF,IP1,IPHES,SCAL,ERR,
1247     &          FIRST,REJECT,FAC1,RPAR,IPAR)
1248C       --- COMPUTE FINITE DIFFERENCES FOR DENSE OUTPUT
1249      IF (ERR.LT.1.D0) THEN
1250         DO I=1,N
1251            Y(I)=Y(I)+ZZ(I+N4)
1252            CONT(I+N5)=ZZ(I)/C(1)
1253         END DO
1254         DO K=1,NS-1
1255            FACT=1.D0/(C(NS-K)-C(NS-K+1))
1256            DO I=1,N
1257               CONT(I+K*N)=(ZZ(I+(NS-K-1)*N)-ZZ(I+(NS-K)*N))*FACT
1258            END DO
1259         END DO
1260         DO J=2,NS
1261            DO K=NS,J,-1
1262               FACT=1.D0/(C(NS-K)-C(NS-K+J))
1263               DO I=1,N
1264                  CONT(I+K*N)=(CONT(I+K*N)-CONT(I+(K-1)*N))*FACT
1265               END DO
1266            END DO
1267         END DO
1268      END IF
1269C *** *** *** *** *** *** ***
1270C *** *** *** *** *** *** ***
1271      ELSE
1272      IF (NS.EQ.7) THEN
1273C *** *** *** *** *** *** ***
1274C *** *** *** *** *** *** ***
1275      IF (FIRST.OR.STARTN.OR.CHANGE) THEN
1276         DO I=1,NNS
1277            ZZ(I)=0.D0
1278            FF(I)=0.D0
1279         END DO
1280      ELSE
1281         HQUOT=H/HOLD
1282         DO K=1,NS
1283            CCQ=C(K)*HQUOT
1284            DO I=1,N
1285               VAL=CONT(I+NS*N)
1286               DO L=NS-1,1,-1
1287                  VAL=CONT(I+L*N)+(CCQ-C(NS-L)+1.D0)*VAL
1288               END DO
1289               ZZ(I+(K-1)*N)=CCQ*VAL
1290            END DO
1291         END DO
1292         DO I=1,N
1293            Z1I=ZZ(I)
1294            Z2I=ZZ(I+N)
1295            Z3I=ZZ(I+N2)
1296            Z4I=ZZ(I+N3)
1297            Z5I=ZZ(I+N4)
1298            Z6I=ZZ(I+N5)
1299            Z7I=ZZ(I+N6)
1300            FF(I)=TI711*Z1I+TI712*Z2I+TI713*Z3I+TI714*Z4I+TI715*Z5I
1301     *               +TI716*Z6I+TI717*Z7I
1302            FF(I+N)=TI721*Z1I+TI722*Z2I+TI723*Z3I+TI724*Z4I+TI725*Z5I
1303     *               +TI726*Z6I+TI727*Z7I
1304            FF(I+N2)=TI731*Z1I+TI732*Z2I+TI733*Z3I+TI734*Z4I+TI735*Z5I
1305     *               +TI736*Z6I+TI737*Z7I
1306            FF(I+N3)=TI741*Z1I+TI742*Z2I+TI743*Z3I+TI744*Z4I+TI745*Z5I
1307     *               +TI746*Z6I+TI747*Z7I
1308            FF(I+N4)=TI751*Z1I+TI752*Z2I+TI753*Z3I+TI754*Z4I+TI755*Z5I
1309     *               +TI756*Z6I+TI757*Z7I
1310            FF(I+N5)=TI761*Z1I+TI762*Z2I+TI763*Z3I+TI764*Z4I+TI765*Z5I
1311     *               +TI766*Z6I+TI767*Z7I
1312            FF(I+N6)=TI771*Z1I+TI772*Z2I+TI773*Z3I+TI774*Z4I+TI775*Z5I
1313     *               +TI776*Z6I+TI777*Z7I
1314         END DO
1315      END IF
1316C *** *** *** *** *** *** ***
1317C  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
1318C *** *** *** *** *** *** ***
1319      NEWT=0
1320      NIT=NIT1+10
1321      EXPMI=1.0D0/EXPMNS
1322      FNEWT=MAX(10*UROUND/RTOL1,MIN(0.03D0,RTOL1**(EXPMI-1.0D0)))
1323      FACCON=MAX(FACCON,UROUND)**0.8D0
1324      THETA=ABS(THET)
1325 240  CONTINUE
1326      IF (NEWT.GE.NIT) GOTO 78
1327C ---     COMPUTE THE RIGHT-HAND SIDE
1328      DO K=0,NS-1
1329         IADD=K*N
1330         DO I=1,N
1331            CONT(I)=Y(I)+ZZ(IADD+I)
1332         END DO
1333         CALL FCN(N,X+C(K+1)*H,CONT,ZZ(IADD+1),RPAR,IPAR)
1334      END DO
1335      NFCN=NFCN+NS
1336C ---     SOLVE THE LINEAR SYSTEMS
1337      DO I=1,N
1338         Z1I=ZZ(I)
1339         Z2I=ZZ(I+N)
1340         Z3I=ZZ(I+N2)
1341         Z4I=ZZ(I+N3)
1342         Z5I=ZZ(I+N4)
1343         Z6I=ZZ(I+N5)
1344         Z7I=ZZ(I+N6)
1345         ZZ(I)=TI711*Z1I+TI712*Z2I+TI713*Z3I+TI714*Z4I+TI715*Z5I
1346     *            +TI716*Z6I+TI717*Z7I
1347         ZZ(I+N)=TI721*Z1I+TI722*Z2I+TI723*Z3I+TI724*Z4I+TI725*Z5I
1348     *            +TI726*Z6I+TI727*Z7I
1349         ZZ(I+N2)=TI731*Z1I+TI732*Z2I+TI733*Z3I+TI734*Z4I+TI735*Z5I
1350     *            +TI736*Z6I+TI737*Z7I
1351         ZZ(I+N3)=TI741*Z1I+TI742*Z2I+TI743*Z3I+TI744*Z4I+TI745*Z5I
1352     *            +TI746*Z6I+TI747*Z7I
1353         ZZ(I+N4)=TI751*Z1I+TI752*Z2I+TI753*Z3I+TI754*Z4I+TI755*Z5I
1354     *            +TI756*Z6I+TI757*Z7I
1355         ZZ(I+N5)=TI761*Z1I+TI762*Z2I+TI763*Z3I+TI764*Z4I+TI765*Z5I
1356     *            +TI766*Z6I+TI767*Z7I
1357         ZZ(I+N6)=TI771*Z1I+TI772*Z2I+TI773*Z3I+TI774*Z4I+TI775*Z5I
1358     *            +TI776*Z6I+TI777*Z7I
1359      END DO
1360      CALL SLVRAR(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1361     &        M1,M2,NM1,FAC1,E1,LDE1,ZZ,FF,IP1,IPHES,IER,IJOB)
1362      DO K=1,3
1363         IAD=(K-1)*2*NM1+1
1364         CALL SLVRAI(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1365     &      M1,M2,NM1,ALPHN(K),BETAN(K),EE2(1,IAD),EE2(1,IAD+NM1),LDE1,
1366     &      ZZ(1+(2*K-1)*N),ZZ(1+2*K*N),FF(1+(2*K-1)*N),
1367     &      FF(1+2*K*N),CONT,IP2((K-1)*NM1+1),IPHES,IER,IJOB)
1368      END DO
1369      NSOL=NSOL+1
1370      NEWT=NEWT+1
1371      DYNO=0.D0
1372      DO I=1,N
1373         DENOM=SCAL(I)
1374         DO K=0,NS-1
1375            DYNO=DYNO+(ZZ(I+K*N)/DENOM)**2
1376         END DO
1377      END DO
1378      DYNO=DSQRT(DYNO/NNS)
1379C ---     BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE
1380      IF (NEWT.GT.1.AND.NEWT.LT.NIT) THEN
1381         THQ=DYNO/DYNOLD
1382         IF (NEWT.EQ.2) THEN
1383            THETA=THQ
1384         ELSE
1385            THETA=SQRT(THQ*THQOLD)
1386         END IF
1387         THQOLD=THQ
1388         IF (THETA.LT.0.99D0) THEN
1389            FACCON=THETA/(1.0D0-THETA)
1390            DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)/FNEWT
1391            IF (DYTH.GE.1.0D0) THEN
1392               QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH))
1393               HHFAC=0.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT))
1394               H=HHFAC*H
1395               REJECT=.TRUE.
1396               LAST=.FALSE.
1397               IF (HHFAC.LE.0.5D0) UNEXN=.TRUE.
1398               IF (CALJAC) GOTO 20
1399               GOTO 10
1400            END IF
1401         ELSE
1402            GOTO 78
1403         END IF
1404      END IF
1405      DYNOLD=MAX(DYNO,UROUND)
1406      DO I=1,N
1407         Z1I=FF(I)+ZZ(I)
1408         Z2I=FF(I+N)+ZZ(I+N)
1409         Z3I=FF(I+N2)+ZZ(I+N2)
1410         Z4I=FF(I+N3)+ZZ(I+N3)
1411         Z5I=FF(I+N4)+ZZ(I+N4)
1412         Z6I=FF(I+N5)+ZZ(I+N5)
1413         Z7I=FF(I+N6)+ZZ(I+N6)
1414         FF(I)=Z1I
1415         FF(I+N)=Z2I
1416         FF(I+N2)=Z3I
1417         FF(I+N3)=Z4I
1418         FF(I+N4)=Z5I
1419         FF(I+N5)=Z6I
1420         FF(I+N6)=Z7I
1421         ZZ(I)=T711*Z1I+T712*Z2I+T713*Z3I+T714*Z4I+T715*Z5I
1422     *            +T716*Z6I+T717*Z7I
1423         ZZ(I+N)=T721*Z1I+T722*Z2I+T723*Z3I+T724*Z4I+T725*Z5I
1424     *            +T726*Z6I+T727*Z7I
1425         ZZ(I+N2)=T731*Z1I+T732*Z2I+T733*Z3I+T734*Z4I+T735*Z5I
1426     *            +T736*Z6I+T737*Z7I
1427         ZZ(I+N3)=T741*Z1I+T742*Z2I+T743*Z3I+T744*Z4I+T745*Z5I
1428     *            +T746*Z6I+T747*Z7I
1429         ZZ(I+N4)=T751*Z1I+T752*Z2I+T753*Z3I+T754*Z4I+T755*Z5I
1430     *            +T756*Z6I+T757*Z7I
1431         ZZ(I+N5)=T761*Z1I+T762*Z2I+T763*Z3I+T764*Z4I+T765*Z5I
1432     *            +T766*Z6I+T767*Z7I
1433         ZZ(I+N6)=T771*Z1I+Z2I+Z4I+Z6I
1434      END DO
1435      IF (FACCON*DYNO.GT.FNEWT) GOTO 240
1436C --- ERROR ESTIMATION
1437      CALL ESTRAV (N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1438     &          H,DD,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,NS,NNS,E1,LDE1,
1439     &          ZZ,CONT,FF,IP1,IPHES,SCAL,ERR,
1440     &          FIRST,REJECT,FAC1,RPAR,IPAR)
1441C       --- COMPUTE FINITE DIFFERENCES FOR DENSE OUTPUT
1442      IF (ERR.LT.1.D0) THEN
1443         DO I=1,N
1444            Y(I)=Y(I)+ZZ(I+(NS-1)*N)
1445            CONT(I+NS*N)=ZZ(I)/C(1)
1446         END DO
1447         DO K=1,NS-1
1448            FACT=1.D0/(C(NS-K)-C(NS-K+1))
1449            DO I=1,N
1450               CONT(I+K*N)=(ZZ(I+(NS-K-1)*N)-ZZ(I+(NS-K)*N))*FACT
1451            END DO
1452         END DO
1453         DO J=2,NS
1454            DO K=NS,J,-1
1455               FACT=1.D0/(C(NS-K)-C(NS-K+J))
1456               DO I=1,N
1457                  CONT(I+K*N)=(CONT(I+K*N)-CONT(I+(K-1)*N))*FACT
1458               END DO
1459            END DO
1460         END DO
1461      END IF
1462C *** *** *** *** *** *** ***
1463C *** *** *** *** *** *** ***
1464      ELSE
1465CASE       (NS.EQ.1)
1466C *** *** *** *** *** *** ***
1467C *** *** *** *** *** *** ***
1468      IF (FIRST.OR.STARTN.OR.CHANGE) THEN
1469         DO I=1,NS
1470            ZZ(I)=0.D0
1471            FF(I)=0.D0
1472         END DO
1473      ELSE
1474         HQUOT=H/HOLD
1475         DO I=1,N
1476            Z1I=HQUOT*CONT(I+N)
1477            ZZ(I)=Z1I
1478            FF(I)=Z1I
1479         END DO
1480      END IF
1481C *** *** *** *** *** *** ***
1482C  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
1483C *** *** *** *** *** *** ***
1484      NEWT=0
1485      NIT=NIT1-3
1486      EXPMI=1.0D0/EXPMNS
1487      FNEWT=MAX(10*UROUND/RTOL1,0.03D0)
1488      FACCON=MAX(FACCON,UROUND)**0.8D0
1489      THETA=ABS(THET)
1490 440  CONTINUE
1491      IF (NEWT.GE.NIT) GOTO 78
1492C ---     COMPUTE THE RIGHT-HAND SIDE
1493      DO I=1,N
1494         CONT(I)=Y(I)+ZZ(I)
1495      END DO
1496      CALL FCN(N,XPH,CONT,ZZ,RPAR,IPAR)
1497      NFCN=NFCN+1
1498C ---     SOLVE THE LINEAR SYSTEMS
1499      CALL SLVRAR(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1500     &          M1,M2,NM1,FAC1,E1,LDE1,ZZ,FF,IP1,IPHES,IER,IJOB)
1501      NSOL=NSOL+1
1502      NEWT=NEWT+1
1503      DYNO=0.D0
1504      DO I=1,N
1505         DENOM=SCAL(I)
1506         DYNO=DYNO+(ZZ(I)/DENOM)**2
1507      END DO
1508      DYNO=DSQRT(DYNO/NNS)
1509C ---     BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE
1510      IF (NEWT.GT.1.AND.NEWT.LT.NIT) THEN
1511         THQ=DYNO/DYNOLD
1512         IF (NEWT.EQ.2) THEN
1513            THETA=THQ
1514         ELSE
1515            THETA=SQRT(THQ*THQOLD)
1516         END IF
1517         THQOLD=THQ
1518         IF (THETA.LT.0.99D0) THEN
1519            FACCON=THETA/(1.0D0-THETA)
1520            DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)/FNEWT
1521            IF (DYTH.GE.1.0D0) THEN
1522               QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH))
1523               HHFAC=0.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT))
1524               H=HHFAC*H
1525               REJECT=.TRUE.
1526               LAST=.FALSE.
1527               IF (HHFAC.LE.0.5D0) UNEXN=.TRUE.
1528               IF (CALJAC) GOTO 20
1529               GOTO 10
1530            END IF
1531         ELSE
1532            GOTO 78
1533         END IF
1534      END IF
1535      DYNOLD=MAX(DYNO,UROUND)
1536      DO I=1,N
1537         F1I=FF(I)+ZZ(I)
1538         FF(I)=F1I
1539         ZZ(I)=F1I
1540      END DO
1541      IF (FACCON*DYNO.GT.FNEWT) GOTO 440
1542C --- ERROR ESTIMATION
1543      CALL ESTRAV (N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1544     &          H,DD,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,NS,NNS,E1,LDE1,
1545     &          ZZ,CONT,FF,IP1,IPHES,SCAL,ERR,
1546     &          FIRST,REJECT,FAC1,RPAR,IPAR)
1547C       --- COMPUTE FINITE DIFFERENCES FOR DENSE OUTPUT
1548      IF (ERR.LT.1.D0) THEN
1549         DO I=1,N
1550            Y(I)=Y(I)+ZZ(I)
1551            CONT(I+N)=ZZ(I)
1552         END DO
1553      END IF
1554C *** *** *** *** *** *** ***
1555C *** *** *** *** *** *** ***
1556      END IF
1557      END IF
1558      END IF
1559C *** *** *** *** *** *** ***
1560C *** *** *** *** *** *** ***
1561C --- COMPUTATION OF HNEW
1562C --- WE REQUIRE .2<=HNEW/H<=8.
1563      FAC=MIN(SAFE,(1+2*NIT)*SAFE/(NEWT+2*NIT))
1564      QUOT=MAX(FACR,MIN(FACL,ERR**EXPO/FAC))
1565      HNEW=H/QUOT
1566C *** *** *** *** *** *** ***
1567C  IS THE ERROR SMALL ENOUGH ?
1568C *** *** *** *** *** *** ***
1569      IF (ERR.LT.1.D0) THEN
1570C --- STEP IS ACCEPTED
1571         FIRST=.FALSE.
1572         NACCPT=NACCPT+1
1573         IF (PRED.AND..NOT.CHANGE) THEN
1574C       --- PREDICTIVE CONTROLLER OF GUSTAFSSON
1575            IF (NACCPT.GT.1) THEN
1576               FACGUS=(HACC/H)*(ERR**2/ERRACC)**EXPO/SAFE
1577               FACGUS=MAX(FACR,MIN(FACL,FACGUS))
1578               QUOT=MAX(QUOT,FACGUS)
1579               HNEW=H/QUOT
1580            END IF
1581            HACC=H
1582            ERRACC=MAX(1.0D-2,ERR)
1583         END IF
1584         XOLD=X
1585         HOLD=H
1586         X=XPH
1587C       --- UPDATE SCALING
1588         IF (ITOL.EQ.0) THEN
1589            DO I=1,N
1590               SCAL(I)=ATOL1+RTOL1*ABS(Y(I))
1591            END DO
1592         ELSE
1593            DO I=1,N
1594               QUOTT=ATOL(I)/RTOL(I)
1595               RTOL1=0.1D0*RTOL(I)**EXPMNS
1596               ATOL1=RTOL1*QUOTT
1597               SCAL(I)=ATOL1+RTOL1*ABS(Y(I))
1598            END DO
1599         END IF
1600         IF (IOUT.NE.0) THEN
1601             NRSOL=NACCPT+1
1602             XSOL=X
1603             XOSOL=XOLD
1604             DO I=1,N
1605                CONT(I)=Y(I)
1606             END DO
1607             NSOLU=N
1608             HSOL=HOLD
1609             CALL SOLOUT(NRSOL,XOSOL,XSOL,Y,CONT,LRC,NSOLU,
1610     &                   RPAR,IPAR,IRTRN)
1611             IF (IRTRN.LT.0) GOTO 179
1612         END IF
1613         CALJAC=.FALSE.
1614         IF (LAST) THEN
1615            H=HOPT
1616            IDID=1
1617            RETURN
1618         END IF
1619         CALL FCN(N,X,Y,Y0,RPAR,IPAR)
1620         NFCN=NFCN+1
1621         HNEW=POSNEG*MIN(ABS(HNEW),HMAXN)
1622         HOPT=HNEW
1623         HOPT=MIN(H,HNEW)
1624         IF (REJECT) HNEW=POSNEG*MIN(ABS(HNEW),ABS(H))
1625         REJECT=.FALSE.
1626         IF ((X+HNEW/QUOT1-XEND)*POSNEG.GE.0.D0) THEN
1627            H=XEND-X
1628            LAST=.TRUE.
1629         ELSE
1630            QT=HNEW/H
1631            HHFAC=H
1632            IF (THETA.LE.THET.AND.QT.GE.QUOT1.AND.QT.LE.QUOT2) THEN
1633               IKEEP=1
1634               GOTO 30
1635            END IF
1636            H=HNEW
1637         END IF
1638         HHFAC=H
1639         IF (THETA.LE.THET) GOTO 20
1640         GOTO 10
1641      ELSE
1642C --- STEP IS REJECTED
1643         REJECT=.TRUE.
1644         LAST=.FALSE.
1645         IF (FIRST) THEN
1646             H=H*0.1D0
1647             HHFAC=0.1D0
1648         ELSE
1649             HHFAC=HNEW/H
1650             H=HNEW
1651         END IF
1652         IF (NACCPT.GE.1) NREJCT=NREJCT+1
1653         IF (CALJAC) GOTO 20
1654         GOTO 10
1655      END IF
1656C --- UNEXPECTED STEP-REJECTION
1657  78  CONTINUE
1658      UNEXP=.TRUE.
1659      IF (IER.NE.0) THEN
1660          NSING=NSING+1
1661          IF (NSING.GE.5) GOTO 176
1662      END IF
1663      H=H*0.5D0
1664      HHFAC=0.5D0
1665      REJECT=.TRUE.
1666      LAST=.FALSE.
1667      IF (CALJAC) GOTO 20
1668      GOTO 10
1669C --- FAIL EXIT
1670 176  CONTINUE
1671      WRITE(6,979)X
1672      WRITE(6,*) ' MATRIX IS REPEATEDLY SINGULAR, IER=',IER
1673      IDID=-4
1674      RETURN
1675 177  CONTINUE
1676      WRITE(6,979)X
1677      WRITE(6,*) ' STEP SIZE T0O SMALL, H=',H
1678      IDID=-3
1679      RETURN
1680 178  CONTINUE
1681      WRITE(6,979)X
1682      WRITE(6,*) ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED'
1683      IDID=-2
1684      RETURN
1685C --- EXIT CAUSED BY SOLOUT
1686 179  CONTINUE
1687      WRITE(6,979)X
1688 979  FORMAT(' EXIT OF RADAU AT X=',E18.4)
1689      IDID=2
1690      RETURN
1691      END
1692C
1693C     END OF SUBROUTINE RADCOV
1694C
1695C ***********************************************************
1696C
1697      SUBROUTINE COERTV(NSMAX)
1698      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1699      COMMON /COE3/ T311,T312,T313,T321,T322,T323,T331,
1700     *    TI311,TI312,TI313,TI321,TI322,TI323,TI331,TI332,TI333
1701      COMMON /COE5/ T511,T512,T513,T514,T515,T521,T522,T523,T524,T525,
1702     *    T531,T532,T533,T534,T535,T541,T542,T543,T544,T545,T551,
1703     *    TI511,TI512,TI513,TI514,TI515,TI521,TI522,TI523,TI524,TI525,
1704     *    TI531,TI532,TI533,TI534,TI535,TI541,TI542,TI543,TI544,TI545,
1705     *    TI551,TI552,TI553,TI554,TI555
1706      COMMON /COE7/ T711,T712,T713,T714,T715,T716,T717,
1707     *              T721,T722,T723,T724,T725,T726,T727,
1708     *              T731,T732,T733,T734,T735,T736,T737,
1709     *              T741,T742,T743,T744,T745,T746,T747,
1710     *              T751,T752,T753,T754,T755,T756,T757,
1711     *              T761,T762,T763,T764,T765,T766,T767,T771,
1712     *              TI711,TI712,TI713,TI714,TI715,TI716,TI717,
1713     *              TI721,TI722,TI723,TI724,TI725,TI726,TI727,
1714     *              TI731,TI732,TI733,TI734,TI735,TI736,TI737,
1715     *              TI741,TI742,TI743,TI744,TI745,TI746,TI747,
1716     *              TI751,TI752,TI753,TI754,TI755,TI756,TI757,
1717     *              TI761,TI762,TI763,TI764,TI765,TI766,TI767,
1718     *              TI771,TI772,TI773,TI774,TI775,TI776,TI777
1719C ---
1720      T311 =  0.9123239487089294279155D-01
1721      T312 = -0.1412552950209542084280D+00
1722      T313 = -0.3002919410514742449186D-01
1723      T321 =  0.2417179327071070189575D+00
1724      T322 =  0.2041293522937999319960D+00
1725      T323 =  0.3829421127572619377954D+00
1726      T331 =  0.9660481826150929361906D+00
1727      TI311 =  0.4325579890063155351024D+01
1728      TI312 =  0.3391992518158098695428D+00
1729      TI313 =  0.5417705399358748711865D+00
1730      TI321 = -0.4178718591551904727346D+01
1731      TI322 = -0.3276828207610623870825D+00
1732      TI323 =  0.4766235545005504519601D+00
1733      TI331 = -0.5028726349457868759512D+00
1734      TI332 =  0.2571926949855605429187D+01
1735      TI333 = -0.5960392048282249249688D+00
1736      IF (NSMAX.LE.3) RETURN
1737      T511 = -0.1251758622050104589014D-01
1738      T512 = -0.1024204781790882707009D-01
1739      T513 =  0.4767387729029572386318D-01
1740      T514 = -0.1147851525522951470794D-01
1741      T515 = -0.1401985889287541028108D-01
1742      T521 = -0.1491670151895382429004D-02
1743      T522 =  0.5017286451737105816299D-01
1744      T523 = -0.9433181918161143698066D-01
1745      T524 = -0.7668830749180162885157D-02
1746      T525 =  0.2470857842651852681253D-01
1747      T531 = -0.7298187638808714862266D-01
1748      T532 = -0.2305395340434179467214D+00
1749      T533 =  0.1027030453801258997922D+00
1750      T534 =  0.1939846399882895091122D-01
1751      T535 =  0.8180035370375117083639D-01
1752      T541 = -0.3800914400035681041264D+00
1753      T542 =  0.3778939022488612495439D+00
1754      T543 =  0.4667441303324943592896D+00
1755      T544 =  0.4076011712801990666217D+00
1756      T545 =  0.1996824278868025259365D+00
1757      T551 = -0.9219789736812104884883D+00
1758      TI511 = -0.3004156772154440162771D+02
1759      TI512 = -0.1386510785627141316518D+02
1760      TI513 = -0.3480002774795185561828D+01
1761      TI514 =  0.1032008797825263422771D+01
1762      TI515 = -0.8043030450739899174753D+00
1763      TI521 =  0.5344186437834911598895D+01
1764      TI522 =  0.4593615567759161004454D+01
1765      TI523 = -0.3036360323459424298646D+01
1766      TI524 =  0.1050660190231458863860D+01
1767      TI525 = -0.2727786118642962705386D+00
1768      TI531 =  0.3748059807439804860051D+01
1769      TI532 = -0.3984965736343884667252D+01
1770      TI533 = -0.1044415641608018792942D+01
1771      TI534 =  0.1184098568137948487231D+01
1772      TI535 = -0.4499177701567803688988D+00
1773      TI541 = -0.3304188021351900000806D+02
1774      TI542 = -0.1737695347906356701945D+02
1775      TI543 = -0.1721290632540055611515D+00
1776      TI544 = -0.9916977798254264258817D-01
1777      TI545 =  0.5312281158383066671849D+00
1778      TI551 = -0.8611443979875291977700D+01
1779      TI552 =  0.9699991409528808231336D+01
1780      TI553 =  0.1914728639696874284851D+01
1781      TI554 =  0.2418692006084940026427D+01
1782      TI555 = -0.1047463487935337418694D+01
1783      IF (NSMAX.LE.5) RETURN
1784      T711 = -0.2153754627310526422828D-02
1785      T712 =  0.2156755135132077338691D-01
1786      T713 =  0.8783567925144144407326D-02
1787      T714 = -0.4055161452331023898198D-02
1788      T715 =  0.4427232753268285479678D-02
1789      T716 = -0.1238646187952874056377D-02
1790      T717 = -0.2760617480543852499548D-02
1791      T721 =  0.1600025077880428526831D-02
1792      T722 = -0.3813164813441154669442D-01
1793      T723 = -0.2152556059400687552385D-01
1794      T724 =  0.8415568276559589237177D-02
1795      T725 = -0.4031949570224549492304D-02
1796      T726 = -0.6666635339396338181761D-04
1797      T727 =  0.3185474825166209848748D-02
1798      T731 = -0.4059107301947683091650D-02
1799      T732 =  0.5739650893938171539757D-01
1800      T733 =  0.5885052920842679105612D-01
1801      T734 = -0.8560431061603432060177D-02
1802      T735 = -0.6923212665023908924141D-02
1803      T736 = -0.2352180982943338340535D-02
1804      T737 =  0.4169077725297562691409D-03
1805      T741 = -0.1575048807937684420346D-01
1806      T742 = -0.3821469359696835048464D-01
1807      T743 = -0.1657368112729438512412D+00
1808      T744 = -0.3737124230238445741907D-01
1809      T745 =  0.8239007298507719404499D-02
1810      T746 =  0.3115071152346175252726D-02
1811      T747 =  0.2511660491343882192836D-01
1812      T751 = -0.1129776610242208076086D+00
1813      T752 = -0.2491742124652636863308D+00
1814      T753 =  0.2735633057986623212132D+00
1815      T754 =  0.5366761379181770094279D-02
1816      T755 =  0.1932111161012620144312D+00
1817      T756 =  0.1017177324817151468081D+00
1818      T757 =  0.9504502035604622821039D-01
1819      T761 = -0.4583810431839315010281D+00
1820      T762 =  0.5315846490836284292051D+00
1821      T763 =  0.4863228366175728940567D+00
1822      T764 =  0.5265742264584492629141D+00
1823      T765 =  0.2755343949896258141929D+00
1824      T766 =  0.5217519452747652852946D+00
1825      T767 =  0.1280719446355438944141D+00
1826      T771 = -0.8813915783538183763135D+00
1827      TI711 = -0.2581319263199822292761D+03
1828      TI712 = -0.1890737630813985089520D+03
1829      TI713 = -0.4908731481793013119445D+02
1830      TI714 = -0.4110647469661428418112D+01
1831      TI715 = -0.4053447889315563304175D+01
1832      TI716 =  0.3112755366607346076554D+01
1833      TI717 = -0.1646774913558444650169D+01
1834      TI721 = -0.3007390169451292131731D+01
1835      TI722 = -0.1101586607876577132911D+02
1836      TI723 =  0.1487799456131656281486D+01
1837      TI724 =  0.2130388159559282459432D+01
1838      TI725 = -0.1816141086817565624822D+01
1839      TI726 =  0.1134325587895161100083D+01
1840      TI727 = -0.4146990459433035319930D+00
1841      TI731 = -0.8441963188321084681757D+01
1842      TI732 = -0.6505252740575150028169D+00
1843      TI733 =  0.6940670730369876478804D+01
1844      TI734 = -0.3205047525597898431565D+01
1845      TI735 =  0.1071280943546478589783D+01
1846      TI736 = -0.3548507491216221879730D+00
1847      TI737 =  0.9198549132786554154409D-01
1848      TI741 =  0.7467833223502269977153D+02
1849      TI742 =  0.8740858897990081640204D+02
1850      TI743 =  0.4024158737379997877014D+01
1851      TI744 = -0.3714806315158364186639D+01
1852      TI745 = -0.3430093985982317350741D+01
1853      TI746 =  0.2696604809765312378853D+01
1854      TI747 = -0.9386927436075461933568D+00
1855      TI751 =  0.5835652885190657724237D+02
1856      TI752 = -0.1006877395780018096325D+02
1857      TI753 = -0.3036638884256667120811D+02
1858      TI754 = -0.1020020865184865985027D+01
1859      TI755 = -0.1124175003784249621267D+00
1860      TI756 =  0.1890640831000377622800D+01
1861      TI757 = -0.9716486393831482282172D+00
1862      TI761 = -0.2991862480282520966786D+03
1863      TI762 = -0.2430407453687447911819D+03
1864      TI763 = -0.4877710407803786921219D+02
1865      TI764 = -0.2038671905741934405280D+01
1866      TI765 =  0.1673560239861084944268D+01
1867      TI766 = -0.1087374032057106164456D+01
1868      TI767 =  0.9019382492960993738427D+00
1869      TI771 = -0.9307650289743530591157D+02
1870      TI772 =  0.2388163105628114427703D+02
1871      TI773 =  0.3927888073081384382710D+02
1872      TI774 =  0.1438891568549108006988D+02
1873      TI775 = -0.3510438399399361221087D+01
1874      TI776 =  0.4863284885566180701215D+01
1875      TI777 = -0.2246482729591239916400D+01
1876      RETURN
1877      END
1878C
1879C     END OF SUBROUTINE COERTV
1880C
1881C ***********************************************************
1882C
1883      SUBROUTINE COERCV(NS,C,DD,U1,ALPH,BETA)
1884      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1885      DIMENSION C(0:NS),DD(NS),ALPH(NS),BETA(NS)
1886      C(0)=0.D0
1887      C(NS)=1.D0
1888      GOTO (1,11,3,11,5,11,7) NS
1889 11   CONTINUE
1890      RETURN
1891  1   CONTINUE
1892      C(1)=1.0D0
1893      U1=1.0D0
1894      DD(1)=-1.0D0
1895      RETURN
1896  3   CONTINUE
1897      SQ6=DSQRT(6.D0)
1898      C(1)=(4.D0-SQ6)/10.D0
1899      C(2)=(4.D0+SQ6)/10.D0
1900      ST9=9.D0**(1.D0/3.D0)
1901      U1=(6.D0+ST9*(ST9-1))/30.D0
1902      ALP=(12.D0-ST9*(ST9-1))/60.D0
1903      BET=ST9*(ST9+1)*DSQRT(3.D0)/60.D0
1904      CNO=ALP**2+BET**2
1905      U1=1.0D0/U1
1906      ALPH(1)=ALP/CNO
1907      BETA(1)=BET/CNO
1908      RETURN
1909  5   CONTINUE
1910      C(1)=  0.5710419611451768219312D-01
1911      C(2)=  0.2768430136381238276800D+00
1912      C(3)=  0.5835904323689168200567D+00
1913      C(4)=  0.8602401356562194478479D+00
1914      DD(1)= -0.2778093394406463730479D+02
1915      DD(2)=  0.3641478498049213152712D+01
1916      DD(3)= -0.1252547721169118720491D+01
1917      DD(4)=  0.5920031671845428725662D+00
1918      DD(5)= -0.2000000000000000000000D+00
1919      U1=  0.6286704751729276645173D+01
1920      ALPH(1)=  0.3655694325463572258243D+01
1921      BETA(1)=  0.6543736899360077294021D+01
1922      ALPH(2)=  0.5700953298671789419170D+01
1923      BETA(2)=  0.3210265600308549888425D+01
1924      RETURN
1925  7   CONTINUE
1926      C(1)=  0.2931642715978489197205D-01
1927      C(2)=  0.1480785996684842918500D+00
1928      C(3)=  0.3369846902811542990971D+00
1929      C(4)=  0.5586715187715501320814D+00
1930      C(5)=  0.7692338620300545009169D+00
1931      C(6)=  0.9269456713197411148519D+00
1932      DD(1)= -0.5437443689412861451458D+02
1933      DD(2)=  0.7000024004259186512041D+01
1934      DD(3)= -0.2355661091987557192256D+01
1935      DD(4)=  0.1132289066106134386384D+01
1936      DD(5)= -0.6468913267673587118673D+00
1937      DD(6)=  0.3875333853753523774248D+00
1938      DD(7)= -0.1428571428571428571429D+00
1939      U1=  0.8936832788405216337302D+01
1940      ALPH(1)=  0.4378693561506806002523D+01
1941      BETA(1)=  0.1016969328379501162732D+02
1942      ALPH(2)=  0.7141055219187640105775D+01
1943      BETA(2)=  0.6623045922639275970621D+01
1944      ALPH(3)=  0.8511834825102945723051D+01
1945      BETA(3)=  0.3281013624325058830036D+01
1946      RETURN
1947      END
1948C
1949C     END OF SUBROUTINE COERCV
1950C
1951C ***********************************************************
1952C
1953      DOUBLE PRECISION FUNCTION CONTRA(I,X,CONT,LRC)
1954C ----------------------------------------------------------
1955C     THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT. IT PROVIDES AN
1956C     APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT X.
1957C     IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR
1958C     THE LAST SUCCESSFULLY COMPUTED STEP (BY RADAU).
1959C ----------------------------------------------------------
1960      PARAMETER (NSDIM=7)
1961      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1962      DIMENSION CONT(LRC)
1963      COMMON/WEIGHT/NN,NS,XSOL,HSOL,C(0:NSDIM)
1964      S=(X-XSOL)/HSOL+1.D0
1965      CONTRA=CONT(I+NS*NN)
1966      DO K=NS-1,0,-1
1967         CONTRA=CONT(I+K*NN)+(S-C(NS-K))*CONTRA
1968      END DO
1969      RETURN
1970      END
1971C
1972C     END OF FUNCTION CONTRA
1973C
1974C ***********************************************************
1975