      PROGRAM LSODET
C...
C...  THE FOLLOWING PROGRAM CALLS AN ODE APPLICATION, DEFINED AS THE
C...  USER-SUPPLIED SUBROUTINES INITAL, DERV AND PRINT, PLUS DATA.
C...  THE SYSTEM OF ODES PROGRAMMED IN SUBROUTINE DERV IS INTEGRATED
C...  BY LSODE (THE BASIC SOLVER IN ODEPACK, A LIBRARY OF INTEGRATORS
C...  DEVELOPED BY ALAN C. HINDMARSH, LAWRENCE LIVERMORE NATIONAL LAB-
C...  ORATORY).
C...
C...  THE MODEL INITIAL CONDITIONS ARE SET IN SUBROUTINE INITAL, AND
C...  THE MODEL DERIVATIVES ARE PROGRAMMED IN SUBROUTINE DERV.  THE
C...  NUMERICAL SOLUTION IS PRINTED AND PLOTTED IN SUBROUTINE PRINT.
C...
C...  SUBROUTINES INITAL, DERV AND PRINT.  THE FOLLOWING CODING IS FOR
C...  250 ORDINARY DIFFERENTIAL EQUATIONS (ODES).  IF MORE ODES ARE TO
C...  BE INTEGRATED, ALL OF THE 250*S SHOULD BE CHANGED TO THE REQUIRED
C...  NUMBER.  DOUBLE PRECISION IS USED FOR THE INTEGRATION OF STIFF
C...  ODES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    Y(250)
     2      /F/    F(250)
C...
C...  COMMON AREA TO PROVIDE THE INPUT/OUTPUT UNIT NUMBERS TO OTHER
C...  SUBROUTINES
      COMMON/IO/       NI,        NO
C...
C...  ABSOLUTE DIMENSIONING OF THE ARRAYS REQUIRED BY LSODE.  ARRAY
C...  WORK IS SIZED BY THE FOLLOWING FORMULA FOR BANDED PROCESSING
C...  OF THE ODE JACOBIAN MATRIX VIA FINITE DIFFERENCES: 22+17*NEQN+
C...  (2*ML+MU)*NEQN.  FOR NEQN = 250, ML = MU = 3, 22+17*250+(2*3+3)*
C...  250 = 6600 (APPROXIMATELY).  IF FEWER THAN 250 ODES ARE TO BE
C...  INTEGRATED, THEN ML AND MU CAN BE LARGER FOR THE GIVEN SIZE OF
C...  WORK.  IF NEQN GT 250 AND/OR ML AND MU GT 3, WORK MUST BE DIMEN-
C...  SIONED TO A LARGER SIZE
      DIMENSION YV(250), RWORK(6600), IWORK(275)
C...
C...  EXTERNAL THE DERIVATIVE AND ODE JACOBIAN MATRIX ROUTINES CALLED BY
C...  LSODE
      EXTERNAL FCN, JAC
C...
C...  ARRAY FOR THE TITLE (FIRST LINE OF DATA), CHARACTERS  END OF RUNS
      CHARACTER TITLE(20)*4, ENDRUN(3)*4
C...
C...  VARIABLE FOR THE TYPE OF ERROR CRITERION
      CHARACTER*3 ABSREL
C...
C...  DEFINE THE CHARACTERS  END OF RUNS
      DATA ENDRUN/'END ','OF R','UNS '/
C...
C...  DEFINE THE INPUT/OUTPUT UNIT NUMBERS
      NI=5
      NO=6
C...
C...  OPEN THE INPUT/OUTPUT FILES
      OPEN(NI,FILE='DATA',  STATUS='OLD')
      OPEN(NO,FILE='OUTPUT',STATUS='NEW')
C...
C...  INITIALIZE THE RUN COUNTER
      NORUN=0
C...
C...  BEGIN A RUN
1     NORUN=NORUN+1
C...
C...  INITIALIZE THE RUN TERMINATION VARIABLE
      NSTOP=0
C...
C...  READ THE FIRST LINE OF DATA
      READ(NI,1000,END=999)(TITLE(I),I=1,20)
C...
C...  TEST FOR  END OF RUNS  IN THE DATA
      DO 2 I=1,3
      IF(TITLE(I).NE.ENDRUN(I))GO TO 3
2     CONTINUE
C...
C...  AN END OF RUNS HAS BEEN READ, SO TERMINATE EXECUTION
999   STOP
C...
C...  READ THE SECOND LINE OF DATA
3     READ(NI,1001,END=999)T0,TF,TP
C...
C...  READ THE THIRD LINE OF DATA
      READ(NI,1002,END=999)NEQN,ERROR
C...
C...  PRINT A DATA SUMMARY
      WRITE(NO,1003)NORUN,(TITLE(I),I=1,20),T0,TF,TP,NEQN,ERROR
C...
C...  INITIALIZE TIME
      T=T0
C...
C...  SET THE INITIAL CONDITIONS
      CALL INITAL
C...
C...  PRINT THE INITIAL CONDITIONS
      CALL PRINT(NI,NO)
C...
C...  SET THE INITIAL CONDITIONS FOR SUBROUTINE LSODE
      DO 5 I=1,NEQN
      YV(I)=Y(I)
5     CONTINUE
C...
C...  SET THE PARAMETERS FOR SUBROUTINE LSODE
      TV=T0
      ITOL=1
      RTOL=ERROR*1.0D-02
      ATOL=ERROR*1.0D-02
      LRW=6600
      LIW=275
      IOPT=0
      ITASK=1
      ISTATE=1
C...
C...  TWO LSODE OPTIONS, AS SELECTED BY THE RUN NUMBER NORUN, ARE
C...  PROGRAMMED HERE. THE TERMINOLOGY IS
C...
C...     IMPLICIT ADAMS METHOD FOR NONSTIFF PROBLEMS
C...
C...     METHOD BASED ON BACKWARD DIFFERENTIATION FORMULAS (BDF) FOR
C...     STIFF PROBLEMS
C...
C...  EACH OF THESE METHODS (NONSTIFF AND STIFF) GENERATES NONLINEAR
C...  ALGEBRAIC EQUATIONS DURING THE NUMERICAL INTEGRATION WHICH ARE
C...  SOLVED BY ONE OF THE FOLLOWING METHODS
C...
C...     FUNCTIONAL ITERATION (A PICARD ITERATION)
C...
C...     CHORD ITERATION, ESSENTIALLY NEWTON*S METHOD, WITH THE ODE
C...     JACOBIAN MATRIX PARTIAL DERIVATIVES APPROXIMATED BY FINITE
C...     DIFFERENCES.  A FULL, BANDED OR DIAGONAL MATRIX MAY BE USED.
C...     THE BANDED AND DIAGONAL FORMS TYPICALLY ARE APPROXIMATIONS TO
C...     THE FULL MATRIX, I.E., CONTAIN A SUBSET OF THE ELEMENTS OF THE
C...     FULL MATRIX.
C...
C...     (1)  NORUN = 1
C...
C...          IMPLICIT ADAMS METHOD, FUNCTIONAL ITERATION
              IF(NORUN.EQ.1)THEN
              METH=1
              MITER=0
C...
C...     (2)  NORUN = 2
C...
C...          IMPLICIT BDF METHOD, CHORD ITERATION WITH AN INTERNALLY
C...          GENERATED (DIFFERENCE QUOTIENT) FULL JACOBIAN
              ELSE IF(NORUN.EQ.2)THEN
              METH=2
              MITER=2
C...
      END IF
C...
C...  THE METHOD FLAG, MF, IS COMPUTED FROM METH AND MITER
      MF=10*METH+MITER
C...
C...  FOR THE NONSTIFF CASE (NORUN = 1), THE INDEPENDENT VARIABLE IS
C...  ADVANCED BY AN ADDITIVE INCREMENT
4        IF(NORUN.EQ.1)THEN
            TOUT=TV+TP
C...
C...  FOR THE STIFF CASE (NORUN = 2), THE INDEPENDENT VARIABLE IS
C...  ADVANCED BY A MULIPLICATIVE FACTOR
         ELSE IF(NORUN.EQ.2)THEN
            TOUT=TV*10.0D0
         END IF
C...
C...  CALL SUBROUTINE LSODE TO COVER ONE PRINT INTERVAL
      CALL LSODE(FCN,NEQN,YV,TV,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,
     1           IOPT,RWORK,LRW,IWORK,LIW,JAC,MF)
C...
C...  PRINT THE SOLUTION
      T=TV
      DO 6 I=1,NEQN
      Y(I)=YV(I)
6     CONTINUE
      CALL PRINT(NI,NO)
C...
C...  TEST FOR AN ERROR CONDITION
      IF(ISTATE.LT.0)THEN
C...
C...     PRINT A MESSAGE INDICATING AN ERROR CONDITION
         WRITE(NO,1004)ISTATE
C...
C...     GO ON TO THE NEXT RUN
         GO TO 1
      END IF
C...
C...  CHECK FOR A RUN TERMINATION
      IF(NSTOP.NE.0)GO TO 1
C...
C...  MONITOR THE PROGRESS OF THE CALCULATION
      WRITE(*,*)TV,TF
C...
C...  CHECK FOR THE END OF THE RUN
      IF(TV.LT.(TF-0.5D0*TP))GO TO 4
C...
C...  THE CURRENT RUN IS COMPLETE, SO PRINT THE COMPUTATIONAL STAT-
C...  ISTICS FOR LSODE AND GO ON TO THE NEXT RUN
      WRITE(NO,8)RWORK(11),IWORK(14),IWORK(11),IWORK(12),IWORK(13)
8     FORMAT(1H ,//,' LSODE COMPUTATIONAL STATISTICS'              ,//,
     1 ' LAST STEP SIZE                       ',              D10.3,//,
     2 ' LAST ORDER OF THE METHOD             ',                I10,//,
     3 ' TOTAL NUMBER OF STEPS TAKEN          ',                I10,//,
     4 ' NUMBER OF FUNCTION EVALUATIONS       ',                I10,//,
     5 ' NUMBER OF JACOBIAN EVALUATIONS       ',                I10,/)
      GO TO 1
C...
C...  *****************************************************************
C...
C...  FORMATS
C...
1000  FORMAT(20A4)
1001  FORMAT(3D10.0)
1002  FORMAT(I5,20X,D10.0)
1003  FORMAT(1H1,
     1 ' RUN NO. - ',I3,2X,20A4,//,
     2 ' INITIAL T - ',D10.3,//,
     3 '   FINAL T - ',D10.3,//,
     4 '   PRINT T - ',D10.3,//,
     5 ' NUMBER OF DIFFERENTIAL EQUATIONS - ',I3,//,
     6 ' MAXIMUM INTEGRATION ERROR - ',D10.3,//,
     7 1H1)
1004  FORMAT(1H ,//,' ISTATE = ',I3,//,
     1 ' INDICATING AN INTEGRATION ERROR, SO THE CURRENT RUN'     ,/,
     2 ' IS TERMINATED.  PLEASE REFER TO THE DOCUMENTATION FOR'   ,/,
     3 ' SUBROUTINE',//,25X,'LSODE',//,
     4 ' FOR AN EXPLANATION OF THESE ERROR INDICATORS'            )
      END
      SUBROUTINE FCN(NEQN,TV,YV,YDOT)
C...
C...  SUBROUTINE FCN IS AN INTERFACE ROUTINE BETWEEN SUBROUTINES LSODE
C...  AND DERV
C...
C...  COMMON AREA
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(1)
     2      /F/      F(1)
C...
C...  VARIABLE DIMENSION THE DEPENDENT AND DERIVATIVE ARRAYS
      DIMENSION YV(NEQN), YDOT(NEQN)
C...
C...  TRANSFER THE INDEPENDENT VARIABLE, DEPENDENT VARIABLE VECTOR
C...  FOR USE IN SUBROUTINE DERV
      T=TV
      DO 1 I=1,NEQN
      Y(I)=YV(I)
1     CONTINUE
C...
C...  EVALUATE THE DERIVATIVE VECTOR
      CALL DERV
C...
C...  TRANSFER THE DERIVATIVE VECTOR FOR USE BY SUBROUTINE LSODE
      DO 2 I=1,NEQN
      YDOT(I)=F(I)
2     CONTINUE
      RETURN
      END
      SUBROUTINE JAC(NEQ,T,Y,ML,MU,PD,NROWPD)
C...
C...  SUBROUTINE JAC IS CALLED ONLY IF AN OPTION OF LSODE IS SELECTED
C...  FOR WHICH THE USER MUST SUPPLY THE ANALYTICAL JACOBIAN MATRIX OF
C...  THE ODE SYSTEM.  THE PROGRAMMING OF THE ANALYTICAL JACOBIAN IN
C...  SUBROUTINE JAC IS DESCRIBED IN DETAIL IN THE DOCUMENTATION COMMENTS
C...  COMMENTS AT THE BEGINNING OF SUBROUTINE LSODE.
C...
      RETURN
      END
      SUBROUTINE INITAL
C...
C...  KROGH'S ODE PROBLEM - SOLUTION OF A STIFF CASE BY LSODE
C...
C...  KROGH'S PROBLEM (1) IS A NONLINEAR ODE PROBLEM WITH AN EXACT
C...  SOLUTION.  THE NUMBER OF ODES AND THE STIFFNESS OF THESE
C...  EQUATIONS CAN BE VARIED BY THE CHOICE OF PARAMETERS.  IN
C...  THIS CASE, N = 20, AND THE STIFFNESS IS SET BY CONSTANTS
C...  READ FROM A DATA FILE.
C...
C...  (1)  GEAR, C. W., NUMERICAL INITIAL VALUE PROBLEMS IN ORDI-
C...       NARY DIFFERENTIAL EQUATIONS, PRENTICE-HALL, ENGLEWOOD
C...       CLIFFS, NJ, 1971, PP 218-220
C...
C...  THE NONLINEAR ODE PROBLEM PROPOSED BY KROGH (1) IS AN IDEAL
C...  TEST PROBLEM FOR ODE SOFTWARE FOR THE FOLLOWING REASONS:
C...
C...     (1)  THE ORDER OF THE PROBLEM (NUMBER OF ODES) CAN BE
C...          SELECTED BY THE USER (IN THE PRESENT CASE, N = 20).
C...
C...     (2)  THE STIFFNESS OF THE ODES CAN BE SELECTED BY THE USER.
C...          THIS POINT REQUIRES SOME CLARIFICATION SINCE THE ODES
C...          ARE NONLINEAR, BUT STIFFNESS IS A LINEAR CONCEPT BASED
C...          ON EIGENVALUES (A STIFF SYSTEM IS GENERALLY DEFINED AS:
C...
C...                 (PROBLEM TIME)*(LARGEST EIGENVALUE)
C...
C...                       IS MUCH GREATER THAN 1)
C...
C...          GENERALLY, THE CONCEPT OF STIFFNESS APPLIED TO NONLINEAR
C...          PROBLEMS MEANS APPLIED TO THE LINEARIZED VERSION OF THE
C...          NONLINEAR EQUATIONS, FOR WHICH EIGENVALUES CAN BE COM-
C...          PUTED.  HOWEVER, SINCE LINEARIZATION IS PERFORMED AROUND
C...          A BASE POINT, THE EIGENVALUES OF THE LINEARIZED SYSTEM
C...          WILL CHANGE WITH THE BASE POINT, I.E., THE LOCAL EIGEN-
C...          VALUES OF A NONLINEAR SYSTEM ARE NOT CONSTANT.
C...
C...          THE NONLINEAR ODES IN THE KROGH PROBLEM HAVE IN THEIR
C...          ANALYTICAL SOLUTION EXPONENTIALS OF THE FORM EXP(B(I)*T)
C...          WHERE T IS THE INDEPENDENT VARIABLE AND THE B(I), I =
C...          1, 2,..., N ARE PARAMETERS SELECTED BY THE USER.  THUS
C...          THE B(I) ARE ANALOGOUS TO EIGENVALUES, BUT THE EXPONENT-
C...          IALS (EXP(B(I)*T)) APPEAR NONLINEARLY IN THE ANALYTICAL
C...          SOLUTION (RATHER THAN AS A LINEAR COMBINATION IN THE
C...          ANALYTICAL SOLUTION OF A SYSTEM OF CONSTANT COEFFICIENT
C...          LINEAR ODES).
C...
C...          IN THIS SENSE, THE STIFFNESS OF THE KROGH TEST PROBLEM
C...          CAN BE SPECIFIED BY THE USER BY SELECTING THE B(I).
C...
C...     (3)  AS THE PRECEDING DISCUSSION INDICATES, THE KROGH TEST
C...          HAS AN ANALYTICAL (EXACT) SOLUTION, WHICH CAN BE USED
C...          AS A STANDARD IN EVALUATING A NUMERICAL INTEGRATOR, I.E.,
C...          THE ERRORS IN NUMERICAL SOLUTIONS CAN BE COMPUTED.
C...
C...  ALL OF THE FEATURES OF THE KROGH TEST PROBLEM ARE ILLUSTRATED IN
C...  THE FOLLOWING CODING FOR TWO CASES:
C...
C...     (1)  NONSTIFF CASE
C...
C...          HERE WE TAKE B(I) = -I, I = 1, 2,..., 20.  THESE CON-
C...          STANTS ARE SET IN DO LOOP 10 IN SUBROUTINE INITAL AT 
C...          THE BEGINNING OF THE NUMERICAL SOLUTION, AND INTEGRATION 
C...          IS BY A NONSTIFF OPTION OF LSODE.
C...
C...     (2)  STIFF CASE
C...
C...          HERE WE TAKE B(I), I = 1, 2,..., 20 TO HAVE THE VALUE
C...          -1000 FOR I LE 10 AND -0.001 FOR I GT 10.  THEY ARE SET 
C...          IN DO LOOP 10 IN SUBROUTINE INITAL AT THE BEGINNING OF THE 
C...          NUMERICAL SOLUTION, AND THE INTEGRATION IS PERFORMED BY A
C...          STIFF OPTION OF LSODE FOR A FULL NUMERICAL JACOBIAN MATRIX.
C...          SINCE THESE VALUES OF B(I) SPAN SIX ORDERS OF MAGNITUDE,
C...          THIS PROBLEM CAN BE CONSIDERED STIFF (AND THEREFORE SHOULD 
C...          NOT BE DONE BY A NONSTIFF (EXPLICIT) INTEGRATOR).
C...
C...  THE MATHEMATICAL DETAILS OF THE PROBLEM ARE GIVEN IN REFERENCE
C...  (1) CITED PREVIOUSLY.  THEY CAN ALSO BE INFERRED FROM THE
C...  FOLLOWING CODING.
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /T/        T,    NSTOP,    NORUN
      COMMON /Y/    Y(20)
      COMMON /F/ DYDT(20)
      COMMON /C/ BETA(20), U(20,20), B(20,20), W(20), N
      COMMON /IO/  NI, NO
C...
C...  READ THE NUMBER OF ODES.  THEN SET THE CONSTANT ARRAY BETA WHICH
C...  DEFINES THE PROBLEM TIME SCALE AND STIFFNESS
      READ(NI,100)N
      DO 10 J=1,N
C...
C...  NONSTIFF CASE, B(I) = -I, I = 1, 2,..., N
      IF(NORUN.EQ.1)THEN
         IF(J.EQ.1)THEN
            BETA(J)=-1.0D0
         ELSE
            BETA(J)=BETA(J-1)-1.0D0
         END IF
C...
C...  STIFF CASE, B(I) = -1000, I LE N/2, B(I) = -0.001, I GT N/2
      ELSE IF(NORUN.EQ.2)THEN
         IF(J.LE.N/2)THEN
            BETA(J)=-1.0D+03
         ELSE
            BETA(J)=-1.0D-03
         END IF
      END IF
  10  CONTINUE
      WRITE(NO,102)N,(BETA(J),J=1,N)
C...
C...  SET THE CONSTANT ARRAY, U
      TERM = 2.0D0/DFLOAT(N)
      DO 20 I=1,N
      DO 30 J=1,N
         U(I,J) = TERM
  30  CONTINUE
         U(I,I) = TERM - 1.0D0
  20  CONTINUE
C...
C...  SET THE CONSTANT ARRAY, B
      DO 40 J=1,N
      DO 50 I=1,N
         SUM = 0.0D0
      DO 60 K=1,N
         SUM = SUM + U(I,K)*BETA(K)*U(K,J)
  60  CONTINUE
         B(I,J) = SUM
  50  CONTINUE
  40  CONTINUE
C...
C...  INITIAL CONDITIONS FROM THE EXACT SOLUTION
      DO 70 J=1,N
         BETAT = BETA(J)*T
C...
C...  AVOID AN OVERFLOW OF THE EXPONENTIAL FUNCTION, EXP(BETAT), IN
C...  THE EXACT SOLUTION
      IF (BETAT.GT.50.0D0) THEN
         W(J) = 0.0D0
C...
C...  AVOID AN UNDERFLOW OF THE EXPONENTIAL FUNCTION, EXP(BETAT), IN
C...  THE EXACT SOLUTION
      ELSE IF (BETAT.LT.-50.0D0) THEN
         W(J) = BETA(J)
      ELSE
         W(J) = BETA(J)/(1.0D0 - (1.0D0 + BETA(J))*DEXP(BETAT))
      END IF
   70 CONTINUE
      DO 80 I=1,N
         SUM = 0.0D0
      DO 90 J=1,N
         SUM = SUM + U(I,J)*W(J)
   90 CONTINUE
         Y(I) = SUM
   80 CONTINUE
C...
C...  PRINT THE OUTPUT HEADING
      WRITE(NO,103)
      WRITE(NO,104)
C...
C...  FORMATS ARE LISTED HERE
 100  FORMAT(I5)
 101  FORMAT(F20.5)
 102  FORMAT(1H ,10X,'N = ',I3,10X,'BETA',20(/,17X,F15.3),//)
 103  FORMAT(1H1)
 104  FORMAT(1H ,6X,'TIME',8X, 'Y(1)',6X, 'YEX(1)',8X, 'PER1',/,
     1                    18X,'Y(20)',5X,'YEX(20)',7X,'PER20',/,
     2                    29X,'PERMAX',/)
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /T/        T,    NSTOP,    NORUN
      COMMON /Y/    Y(20)
      COMMON /F/ DYDT(20)
      COMMON /C/ BETA(20), U(20,20), B(20,20), W(20), N
C...
C...  CALCULATE THE ARRAY, W
      DO 10 I=1,N
         SUM = 0.0D0
      DO 20 J=1,N
         SUM = SUM + U(I,J)*Y(J)
  20  CONTINUE
         W(I) = SUM
  10  CONTINUE
C...
C...  CALCULATE THE N ODE DERIVATIVES
      DO 30 I=1,N
         SUM = 0.0D0
      DO 40 J=1,N
         SUM = SUM + U(I,J)*W(J)*W(J) - B(I,J)*Y(J)
  40  CONTINUE
         DYDT(I) = SUM
  30  CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /T/        T,    NSTOP,    NORUN
      COMMON /Y/    Y(20)
      COMMON /F/ DYDT(20)
      COMMON /C/ BETA(20), U(20,20), B(20,20), W(20), N
      DIMENSION YEX(20)
C...
C...  ABSOLUTE DIMENSIONING OF THE ARRAYS REQUIRED BY SUBROUTINE JMAP
C...  (A, X, XOLD, F, FOLD) AND SUBROUTINE EIGEN (WI, WR, Z, RW, IW).
C...  NOTE THAT ARRAYS A, WR, WI, Z AND RW ARE SINGLE PRECISION IN
C...  ORDER TO CONFORM TO THE CODING IN SUBROUTINES JMAP AND EIGEN
      REAL       A(20,20),   WR(20), WI(20),  Z(20,20), RW(20)
      DIMENSION     X(20), XOLD(20),  F(20),  FOLD(20),PER(20)
      INTEGER      IW(20)
C...
C...  THE FOLLOWING EQUIVALENCE FACILITATES TRANSFER OF THE DEPENDENT
C...  VARIABLE VECTOR IN /Y/ TO ARRAY Y, AND THE DERIVATIVE VECTOR IN
C...  IN /F/ TO ARRAY F FOR USE IN THE SUBSEQUENT CALL TO SUBROUTINE
C...  JMAP
      EQUIVALENCE (Y(1),X(1)),(DYDT(1),F(1))
C...
C...  CALCULATE THE N EXACT SOLUTIONS
      DO 10 J=1,N
         BETAT = BETA(J)*T
C...
C...  AVOID AN OVERFLOW OF THE EXPONENTIAL FUNCTION, EXP(BETAT), IN
C...  THE EXACT SOLUTION
      IF (BETAT.GT.50.0D0) THEN
         W(J) = 0.0D0
C...
C...  AVOID AN UNDERFLOW OF THE EXPONENTIAL FUNCTION, EXP(BETAT), IN
C...  THE EXACT ASOLUTION
      ELSE IF (BETAT.LT.-50.0D0) THEN
         W(J) = BETA(J)
      ELSE
         W(J) = BETA(J)/(1.0D0 - (1.0D0 + BETA(J))*DEXP(BETAT))
      END IF
  10  CONTINUE
      DO 20 I=1,N
         SUM = 0.0D0
      DO 30 J=1,N
         SUM = SUM + U(I,J)*W(J)
  30  CONTINUE
         YEX(I) = SUM
  20  CONTINUE
C...
C...  CALCULATE THE MAXIMUM ERROR
      ERRMAX = DABS((Y(1) - YEX(1))/YEX(1))*100.0D0
      DO 40 J=2,N
         ERR = DABS((Y(J) - YEX(J))/YEX(J))*100.0D0
         IF (ERR.GT.ERRMAX) ERRMAX = ERR
  40  CONTINUE
C...
C...  PRINT THE NUMERICAL AND EXACT SOLUTIONS FOR I = 1 AND N, AND
C...  THE CORRESPONDING ERRORS, AND THE MAXIMUM ERROR
      PER1=DABS((Y(1)-YEX(1))/YEX(1))*100.0D0
      PERN=DABS((Y(N)-YEX(N))/YEX(N))*100.0D0
      PERMAX=ERRMAX
      WRITE(NO,101)T,Y(1),YEX(1),PER1,Y(N),YEX(N),PERN,PERMAX
 101  FORMAT(F11.6,2E12.4,F12.3,/,
     +         11X,2E12.4,F12.3,/,
     +         23X,       F12.3,/)
C...
C...  MAP THE JACOBIAN MATRIX OF THE ODE SYSTEM DEFINED IN SUBROUTINE
C...  DERV, AND COMPUTE ITS TEMPORAL EIGENVALUES
C...  CALL JMAP(N,A,X,XOLD,F,FOLD)
C...
C...  SUBROUTINE EIGEN CALLS A SERIES OF EISPACK ROUTINES TO
C...  COMPUTE THE TEMPORAL EIGENVALUES, AND OPTIONALLY, THE
C...  EIGENVECTORS OF THE 20-ODE SYSTEM JACOBIAN MATRIX
C...  CALL EIGEN(N,A,WR,WI,Z,RW,IW)
      RETURN
      END
