DIMENSION OF THE SYSTEM

NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE VALUE OF F(X,Y): SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR) DOUBLE PRECISION X,Y(N),F(N) F(1)=... ETC. RPAR, IPAR (SEE BELOW)

INITIAL X-VALUE

FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)

INITIAL STEP SIZE GUESS; FOR STIFF EQUATIONS WITH INITIAL TRANSIENT, H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD. THIS CHOICE IS NOT VERY IMPORTANT, THE STEP SIZE IS QUICKLY ADAPTED. (IF H=0.D0, THE CODE PUTS H=1.D-6).

SWITCH FOR RTOL AND ATOL: ITOL=0: BOTH RTOL AND ATOL ARE SCALARS. THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF Y(I) BELOW RTOL*ABS(Y(I))+ATOL ITOL=1: BOTH RTOL AND ATOL ARE VECTORS. THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW RTOL(I)*ABS(Y(I))+ATOL(I).

NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY A DUMMY SUBROUTINE IN THE CASE IJAC=0). FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM SUBROUTINE JAC(N,X,Y,DFY,LDFY,RPAR,IPAR) DOUBLE PRECISION X,Y(N),DFY(LDFY,N) DFY(1,1)= ... LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS FURNISHED BY THE CALLING PROGRAM. IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO BE FULL AND THE PARTIAL DERIVATIVES ARE STORED IN DFY AS DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J) ELSE, THE JACOBIAN IS TAKEN AS BANDED AND THE PARTIAL DERIVATIVES ARE STORED DIAGONAL-WISE AS DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).

SWITCH FOR THE COMPUTATION OF THE JACOBIAN: IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED. IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.

SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN: MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. 0.LE.MLJAC.LT.N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN MATRIX (.GE. NUMBER OF NON-ZERO DIAGONALS BELOW THE MAIN DIAGONAL).

UPPER BANDWITH OF JACOBIAN MATRIX (.GE. NUMBER OF NON- ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). NEED NOT BE DEFINED IF MLJAC=N.

NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS- MATRIX M. IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY MATRIX AND NEEDS NOT TO BE DEFINED; SUPPLY A DUMMY SUBROUTINE IN THIS CASE. IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM SUBROUTINE MAS(N,AM,LMAS,RPAR,IPAR) DOUBLE PRECISION AM(LMAS,N) AM(1,1)= .... IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED AS FULL MATRIX LIKE AM(I,J) = M(I,J) ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED DIAGONAL-WISE AS AM(I-J+MUMAS+1,J) = M(I,J).

GIVES INFORMATION ON THE MASS-MATRIX: IMAS=0: M IS SUPPOSED TO BE THE IDENTITY MATRIX, MAS IS NEVER CALLED. IMAS=1: MASS-MATRIX IS SUPPLIED.

SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX: MLMAS=N: THE FULL MATRIX CASE. THE LINEAR ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. 0.LE.MLMAS.LT.N: MLMAS IS THE LOWER BANDWITH OF THE MATRIX (.GE. NUMBER OF NON-ZERO DIAGONALS BELOW THE MAIN DIAGONAL). MLMAS IS SUPPOSED TO BE .LE. MLJAC.

UPPER BANDWITH OF MASS-MATRIX (.GE. NUMBER OF NON- ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). NEED NOT BE DEFINED IF MLMAS=N. MUMAS IS SUPPOSED TO BE .LE. MUJAC.

NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE NUMERICAL SOLUTION DURING INTEGRATION. IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. SUPPLY A DUMMY SUBROUTINE IF IOUT=0. IT MUST HAVE THE FORM SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N, RPAR,IPAR,IRTRN) DOUBLE PRECISION X,Y(N),CONT(LRC) .... SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH GRID-POINT "X" (THEREBY THE INITIAL VALUE IS THE FIRST GRID-POINT). "XOLD" IS THE PRECEEDING GRID-POINT. "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN IS SET .LT.0, RADAU5 RETURNS TO THE CALLING PROGRAM. ----- CONTINUOUS OUTPUT: ----- DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH THE FUNCTION .GT..GT..GT. CONTR5(I,S,CONT,LRC) .LT..LT..LT. WHICH PROVIDES AN APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE S SHOULD LIE IN THE INTERVAL [XOLD,X]. DO NOT CHANGE THE ENTRIES OF CONT(LRC), IF THE DENSE OUTPUT FUNCTION IS USED.

SWITCH FOR CALLING THE SUBROUTINE SOLOUT: IOUT=0: SUBROUTINE IS NEVER CALLED IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT.

ARRAY OF WORKING SPACE OF LENGTH "LWORK". WORK(1), WORK(2),.., WORK(20) SERVE AS PARAMETERS FOR THE CODE. FOR STANDARD USE OF THE CODE WORK(1),..,WORK(20) MUST BE SET TO ZERO BEFORE CALLING. SEE BELOW FOR A MORE SOPHISTICATED USE. WORK(21),..,WORK(LWORK) SERVE AS WORKING SPACE FOR ALL VECTORS AND MATRICES. "LWORK" MUST BE AT LEAST N*(LJAC+LMAS+3*LE+12)+20 WHERE LJAC=N IF MLJAC=N (FULL JACOBIAN) LJAC=MLJAC+MUJAC+1 IF MLJAC.LT.N (BANDED JAC.) AND LMAS=0 IF IMAS=0 LMAS=N IF IMAS=1 AND MLMAS=N (FULL) LMAS=MLMAS+MUMAS+1 IF MLMAS.LT.N (BANDED MASS-M.) AND LE=N IF MLJAC=N (FULL JACOBIAN) LE=2*MLJAC+MUJAC+1 IF MLJAC.LT.N (BANDED JAC.) IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM STORAGE REQUIREMENT IS LWORK = 4*N*N+12*N+20. IF IWORK(9)=M1.GT.0 THEN "LWORK" MUST BE AT LEAST N*(LJAC+12)+(N-M1)*(LMAS+3*LE)+20 WHERE IN THE DEFINITIONS OF LJAC, LMAS AND LE THE NUMBER N CAN BE REPLACED BY N-M1.

DECLARED LENGTH OF ARRAY "WORK".

INTEGER WORKING SPACE OF LENGTH "LIWORK". IWORK(1),IWORK(2),...,IWORK(20) SERVE AS PARAMETERS FOR THE CODE. FOR STANDARD USE, SET IWORK(1),.., IWORK(20) TO ZERO BEFORE CALLING. IWORK(21),...,IWORK(LIWORK) SERVE AS WORKING AREA. "LIWORK" MUST BE AT LEAST 3*N+20.

DECLARED LENGTH OF ARRAY "IWORK".

REPORTS ON SUCCESSFULNESS UPON RETURN: IDID= 1 COMPUTATION SUCCESSFUL, IDID= 2 COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT) IDID=-1 INPUT IS NOT CONSISTENT, IDID=-2 LARGER NMAX IS NEEDED, IDID=-3 STEP SIZE BECOMES TOO SMALL, IDID=-4 MATRIX IS REPEATEDLY SINGULAR.