      PROGRAM LSODET
C...
C...  AN INTERFACE BETWEEN DSS/2 APPLICATIONS AND INTEGRATOR LSODE
C...
C...  THE FOLLOWING PROGRAM CALLS A DSS/2 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...  TYPE ALL REAL VARIABLES AS DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
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
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    Y(250)
     2      /F/    F(250)
C...
C...  COMMON STAT IS USED TO PASS THE LSODE COMPUTATIONAL STATISTICS 
C...  TO SUBROUTINE PRINT
      COMMON/STAT/HLAST,NST,NFE,NJE,NQU
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 RU','NS  '/
C...
C...  DEFINE THE INPUT/OUTPUT UNIT NUMBERS AND FILES
      NI=5
      NO=6
      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,NMAX,INT,IERROR,ABSREL,ERROR
C...
C...  PRINT A DATA SUMMARY
      WRITE(NO,1003)NORUN,(TITLE(I),I=1,20),
     1                     T0,TF,TP,
     2                     NEQN,NMAX,INT,IERROR,ABSREL,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
      ATOL=1.0D-10
      LRW=6600
      LIW=275
      IOPT=0
      ITASK=1
      ISTATE=1
C...
C...  SEVERAL 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...  SINCE A STIFF ODE PROBLEM IS BEING INTEGRATED, THE STIFF OPTION
C...  WITH INTERNALLY GENERATED JACOBIAN MATRIX IS SELECTED
      NORUN=6
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 ADAMS METHOD, CHORD ITERATION WITH AN INTERNALLY
C...          GENERATED (DIFFERENCE QUOTIENT) FULL JACOBIAN
              ELSE IF(NORUN.EQ.2)THEN
              METH=1
              MITER=2
C...
C...     (3)  NORUN = 3
C...
C...          IMPLICIT ADAMS METHOD, CHORD ITERATION WITH AN INTERNALLY
C...          GENERATED (DIFFERENCE QUOTIENT) DIAGONAL JACOBIAN APPROXI-
C...          TION
              ELSE IF(NORUN.EQ.3)THEN
              METH=1
              MITER=3
C...
C...     (4)  NORUN = 4
C...
C...          IMPLICIT ADAMS METHOD, CHORD ITERATION WITH AN INTERNALLY
C...          GENERATED (DIFFERENCE QUOTIENT) BANDED JACOBIAN
              ELSE IF(NORUN.EQ.4)THEN
              METH=1
              MITER=5
              IWORK(1)=1
              IWORK(2)=1
C...
C...     (5)  NORUN = 5
C...
C...          IMPLICIT BDF METHOD, FUNCTIONAL ITERATION
              ELSE IF(NORUN.EQ.5)THEN
              METH=2
              MITER=0
C...
C...     (6)  NORUN = 6
C...
C...          IMPLICIT BDF METHOD, CHORD ITERATION WITH AN INTERNALLY
C...          GENERATED (DIFFERENCE QUOTIENT) FULL JACOBIAN
              ELSE IF(NORUN.EQ.6)THEN
              METH=2
              MITER=2
C...
C...     (7)  NORUN = 7
C...
C...          IMPLICIT BDF METHOD, CHORD ITERATION WITH AN INTERNALLY
C...          GENERATED (DIFFERENCE QUOTIENT) DIAGONAL JACOBIAN APPROXI-
C...          TION
              ELSE IF(NORUN.EQ.7)THEN
              METH=2
              MITER=3
C...
C...     (8)  NORUN = 8
C...
C...          IMPLICIT BDF METHOD, CHORD ITERATION WITH AN INTERNALLY
C...          GENERATED (DIFFERENCE QUOTIENT) BANDED JACOBIAN
              ELSE IF(NORUN.EQ.8)THEN
              METH=2
              MITER=5
              IWORK(1)=1
              IWORK(2)=1
      END IF
C...
C...  THE METHOD FLAG, MF, IS COMPUTED FROM METH AND MITER
      MF=10*METH+MITER
C...
C...  INITIATE THE INTEGRATION
      TOUT=TV+TP
C...
C...  CALL SUBROUTINE LSODE TO COVER ONE PRINT INTERVAL
4     CALL LSODE(FCN,NEQN,YV,TV,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,
     1           IOPT,RWORK,LRW,IWORK,LIW,JAC,MF)
C...
C...  PASS THE COMPUTATIONAL STATISTICS TO PRINT (VIA COMMON/STAT/)
      HLAST=RWORK(11)
      NST=IWORK(11)
      NFE=IWORK(12)
      NJE=IWORK(13)
      NQU=IWORK(14)
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...  THE INTEGRATION HAS PROCEEDED SATISFACTORILY, SO PREPARE FOR THE
C...  NEXT INTERVAL IN T.  NOTE THAT THE FOLLOWING CALCULATION OF THE
C...  OUTPUT TIME, TOUT, PRODUCES THREE OUTPUT POINTS FOR EACH DECADE
C...  IN T
      TOUT=TV*(1.0D+01)**(1.0D+00/3.0D+00)
C...
C...  CHECK FOR THE END OF THE RUN
      IF(TV.LT.(TF-0.5D+00*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 ,//,' 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(4I5,2X,A3,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 ' PRINT INTERVAL/MINIMUM INTEGRATION INTERVAL - ',I5,//,
     7 ' INTEGRATION ALGORITHM - ',I2,' - LSODE ',//,
     8 ' INTEGRATION ERROR MESSAGES - ',I1,//,
     9 ' ERROR CRITERION - ',A3,//,
     A ' MAXIMUM INTEGRATION ERROR - ',D10.3,//,
     B 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...  DSS/2 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
C...
C...  SUBROUTINE JAC IS A DUMMY ROUTINE TO SATISFY THE LOADER (SINCE JAC
C...  IS DECLARED AS AN EXTERNAL IN THE MAIN PROGRAM LSODET).  JAC IS NOT
C...  ACTUALLY CALLED UNLESS AN OPTION OF LSODE IS SELECTED WHICH REQUIRES
C...  THE USER TO PROVIDE THE ODE ANALYTICAL JACOBIAN MATRIX.  THIS IS
C...  USUALLY NOT THE CASE SINCE FOR MOST PROBLEMS, THE JACOBIAN MATRIX
C...  IS CALCULATED INTERNALLY BY LSODE USING FINITE DIFFERENCES
      RETURN
      END
      SUBROUTINE INITAL
C...
C...  TWO LINEAR CONSTANT COEFFICIENT ORDINARY DIFFERENTIAL EQUATIONS
C...  WITH A STIFFNESS RATIO OF ONE MILLION - PLOTTING OF THE COMPU-
C...  TATIONAL STATISTICS
C...
C...  THE TWO SIMULTANEOUS, LINEAR, CONSTANT COEFFICIENT ORDINARY
C...  DIFFERENTIAL EQUATIONS
C...
C...     DY1/DT = A11*Y1 + A12*Y2                                    (1)
C...
C...     DY2/DT = A21*Y1 + A22*Y2                                    (2)
C...
C...  ARE INTEGRATED SUBJECT TO THE INITIAL CONDITIONS
C...
C...     Y1(0) = 0, Y2(0) = 2                                     (3)(4)
C...
C...  THE ANALYTICAL SOLUTION TO EQUATIONS (1) TO (4) CAN EASILY BE
C...  DERIVED.  FOR EXAMPLE, IF Y1(T) AND Y2(T) ARE ASSUMED TO BE
C...
C...     Y1 = C1*EXP(L*T), Y2 = C2*EXP(L*T)
C...
C...  SUBSTITUTION IN EQUATIONS (1) AND (2) GIVES
C...
C...     L*C1*EXP(L*T) = A11*C1*EXP(L*T) + A12*C2*EXP(L*T)           (5)
C...
C...     L*C2*EXP(L*T) = A21*C1*EXP(L*T) + A22*C2*EX(L*T)            (6)
C...
C...  THE EXPONENTIAL EXP(L*T) CANCELS FROM EQUATIONS (5) AND (6), AND
C...  THE FOLLOWING HOMOGENEOUS LINEAR ALGEBRAIC EQUATIONS FOR C1 AND C2
C...  RESULT
C...
C...     (A11 - L)*C1 +       A12*C2 = 0                             (7)
C...
C...           A21*C1 + (A22 - L)*C2 = 0                             (8)
C...
C...  EQUATIONS (7) AND (8) CONSTITUTE THE CLASSICAL, SECOND-ORDER,
C...  LINEAR ALGERAIC EIGENVALUE PROBLEM WHICH CAN BE SOLVED FOR THE
C...  VALUES OF L THAT DO NOT GIVE THE TRIVIAL SOLUTION C1 = 0, C2 = 0.
C...  THIS WILL BE DONE BY SETTING THE DETERMINANT OF THE COEFFICIENT
C...  MATRIX EQUAL TO ZERO (OTHER ALGORITHMS FOR CALCULATING THE EIGEN-
C...  VALUES, L1 AND L2, ARE AVAILABLE)
C...
C...     (A11 - L)*(A22 - L) - A12*A21 = 0                           (9)
C...
C...  OR
C...
C...     L**2 - (A11 + A22)*L + (A11*A22 - A12*A21) = 0             (10)
C...
C...  SOLUTION BY THE QUADRATIC EQUATION GIVES THE TWO EIGENVALUES, L1
C...  AND L2.  SINCE EQUATIONS (1) AND (2) ARE LINEAR, THE GENERAL
C...  SOLUTION CONSISTS OF A LINEAR COMBINATION OF THE TWO EXPONENTIALS
C...  EXP(L1*T) AND EXP(L2*T), I.E.,
C...
C...     Y1 = C11*EXP(L1*T) + C12*EXP(L2*T)                         (11)
C...
C...     Y2 = C21*EXP(L1*T) + C22*EXP(L2*T)                         (12)
C...
C...  WHERE C11, C12, C21 AND C22 ARE CONSTANTS TO BE DETERMINED (NOTE
C...  THAT (C11,C21) IS THE EIGENVECTOR FOR EIGENVALUE L1, AND (C12,C22)
C...  IS THE EIGENVECTOR FOR L2).  APPLYING INITIAL CONDITIONS (3) AND
C...  (4), PLUS EITHER EQUATION (7) OR (8) FOR L1 AND L2, GIVES THE
C...  RESULT
C...
C...     Y1 = -EXP(L1*T) + EXP(L2*T)                                (13)
C...
C...     Y2 =  EXP(L1*T) + EXP(L2*T)                                (14)
C...
C...  FOR THE SPECIAL CASE A11 = A22 = -A, A12 = A21 = B, L1 = -(A + B),
C...  L2 = -(A - B).  THE FOLLOWING PARTICULAR VALUES OF A AND B ILLUS-
C...  TRATE VARYING DEGREES OF STIFFNESS
C...
C...     CASE 1   A = 50.5        L1 = -100        L1/L2 = 100
C...              B = 49.5        L2 =   -1        (NONSTIFF)
C...
C...     CASE 2   A = 500.5       L1 = -1000       L1/L2 = 1000
C...              B = 499.5       L2 =    -1    (MODERATELY STIFF)
C...
C...     CASE 3   A = 500000.5    L1 = -1000000    L1/L2 = 1000000
C...              B = 499999.5    L2 =       -1     (VERY STIFF)
C...
C...  THE THIRD CASE PROVIDES A STRINGENT TEST OF A NUMERICAL INTEGRATOR
C...  FOR STIFF DIFFERENTIAL EQUATIONS.  FOR EXAMPLE, IF AN EXPLICIT
C...  INTEGRATOR IS USED IN CASE THREE, THE MAXIMUM STEP SIZE DETERMINED
C...  BY THE STABILITY OF THE CALCULATION WILL BE GIVEN APPROXIMATELY BY
C...
C...      ABS(HMAX*L1) LT 2
C...
C...  AND WITH L1 = -1000000, HMAX LT 2/1000000 = 0.000002.  SINCE THE
C...  TIME FOR THE PROBLEM IS DETERMINED BY L2, THROUGH THE EXPONENTIAL
C...  EXP(L2*T) IN EQUATIONS (13) AND (14), THE NUMERICAL SOLUTION
C...  MUST BE COMPUTED APPROXIMATELY OVER THE INTERVAL O LE T LE 10.
C...  THUS, TO COMPUTE THE SOLUTION TO T = 10, 10/0.000002 = 5,000,000
C...  STEPS IN T WILL BE REQUIRED BY THE EXPLICIT INTEGRATOR.  AN
C...  IMPLICIT INTEGRATOR, HOWEVER, WHICH IS NOT LIMITED BY STABILITY
C...  SHOULD PRODUCE A SOLUTION OVER THE INTERVAL 0 LE T LE 10 WITH AT
C...  MOST A FEW HUNDRED STEPS (THE NUMBER OF STEPS  WILL BE DETERMINED
C...  BY THE REQUIRED ACCURACY AND NOT BY THE STABILITY OF THE CALCU-
C...  LATION).  THIS VASTLY IMPROVED EFFICIENCY OF IMPLICIT INTEGRATORS
C...  FOR THIS STIFF PROBLEM IS DEMONSTRATED BY THE FOLLOWING CODING
C...
C...  BECAUSE OF THE STIFFNESS OF THIS PROBLEM, DOUBLE PRECISION CAL-
C...  CULATIONS ARE ADVISED
      DOUBLE PRECISION
     1              T,
     2             Y1,    Y2,
     3          DY1DT, DY2DT,
     4              A,     B
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/    Y1,    Y2
     2      /F/ DY1DT, DY2DT
     3      /C/     A,     B
C...
C...  INITIAL CONDITIONS, EQUATIONS (3) AND (4)
      Y1=0.D+00
      Y2=2.D+00
C...
C...  THE FOLLOWING CONSTANTS IN THE ORDINARY DIFFERENTIAL EQUATIONS
C...  (ODES) CORRESPOND TO A STIFFNESS RATIO OF ONE MILLION (CASE 3
C...  ABOVE)
      A=500000.5D+00
      B=499999.5D+00
      RETURN
      END
      SUBROUTINE DERV
      DOUBLE PRECISION
     1              T,
     2             Y1,    Y2,
     3          DY1DT, DY2DT,
     4              A,     B
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/    Y1,    Y2
     2      /F/ DY1DT, DY2DT
     3      /C/     A,     B
C...
C...  ODES, EQUATIONS (1) AND (2)
      DY1DT = -A*Y1 + B*Y2
      DY2DT =  B*Y1 - A*Y2
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      DOUBLE PRECISION
     1              T,
     2             Y1,    Y2,
     3          DY1DT, DY2DT,
     4              A,     B, HLAST
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/    Y1,    Y2
     2      /F/ DY1DT, DY2DT
     3      /C/     A,     B
C...
C...  COMMON STAT IS USED TO PASS THE LSODE COMPUTATIONAL STATISTICS 
C...  TO SUBROUTINE PRINT
      COMMON/STAT/HLAST,NST,NFE,NJE,NQU
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      REAL         TP(31),  YP(2,31),
     1             P1(31),    P2(31),    P3(31),    P4(31),    P5(31),
     2             P6(31),    P7(31),    P8(31)
C...
C...  INITIALIZE A COUNTER FOR THE PLOTTED STATISTICS
      DATA IP/0/
C...
C...  PRINT A HEADING FOR THE COMPUTATIONAL STATISTICS
      IF(T.EQ.0.D+00)WRITE(NO,2)
2     FORMAT(1H ,//,9X,1HT,
     1 10X,2HHU,9X,3HNST,9X,3HNFE,8X,3HNJE,6X,3HNQU)
C...
C...  INITIALIZE THREE VARIABLES FOR COMPUTING THE NET NUMBER OF STEPS,
C...  FUNCTION EVALUATIONS AND JACOBIAN MATRIX EVALUATIONS
      IF(T.GT.0.D+00)GO TO 3
      NET1=0
      NET2=0
      NET3=0
      RETURN
C...
C...  PRINT THE COMPUTATIONAL STATISTICS
3     NET1=NST-NET1
      NET2=NFE-NET2
      NET3=NJE-NET3
      WRITE(NO,1)T,HU,NST,NET1,NFE,NET2,NJE,NET3,NQU
1     FORMAT(D10.3,D12.4,
     1       I6,2H (,I3,1H),I6,2H (,I3,1H),I5,3H  (,I2,1H),I9)
C...
C...  STORE THE SOLUTION AND COMPUTATIONAL STATISTICS FOR SUBSEQUENT
C...  PLOTTING
      IP=IP+1
      TP(IP)=ALOG10(SNGL(T))
      YP(1,IP)=SNGL(Y1)
      YP(2,IP)=SNGL(Y2)
      P1(IP)=SNGL(HLAST)
      P2(IP)=FLOAT(NST)
      P3(IP)=FLOAT(NET1)
      P4(IP)=FLOAT(NFE)
      P5(IP)=FLOAT(NET2)
      P6(IP)=FLOAT(NJE)
      P7(IP)=FLOAT(NET3)
      P8(IP)=FLOAT(NQU)
      NET1=NST
      NET2=NFE
      NET3=NJE
C...
C...  PLOT THE SOLUTION AND COMPUTATIONAL STATISTICS
      IF(IP.LT.31)RETURN
C...  ******************************************************************
C...
C...  THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE-DEPENDENT
C...  AND WILL HAVE TO BE CHANGED FOR THE USER*S LOCAL COMPUTER
C...
C...  CYBER PLOTTING AT LEHIGH UNIVERSITY
C...
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(5.,-9.,2.,5.,0.,0.4)
C...  CALL QIKPLT(TP,YP,IP,8H*LOG(T)*,14H*Y1(T), Y2(T)*,2H**,2)
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(5.,-9.,2.,5.,0.,-0.0)
C...  CALL QIKPLT(TP,P1,IP,8H*LOG(T)*,12H*INTEG STEP*,2H**,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(5.,-9.,2.,5.,0.,-0.0)
C...  CALL QIKPLT(TP,P2,IP,8H*LOG(T)*,14H*CUM NO STEPS*,2H**,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(5.,-9.,2.,5.,0.,-0.0)
C...  CALL QIKPLT(TP,P3,IP,8H*LOG(T)*,14H*INC NO STEPS*,2H**,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(5.,-9.,2.,5.,0.,-0.0)
C...  CALL QIKPLT(TP,P4,IP,8H*LOG(T)*,14H*CUM NO DERVS*,2H**,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(5.,-9.,2.,5.,0.,-0.0)
C...  CALL QIKPLT(TP,P5,IP,8H*LOG(T)*,14H*INC NO DERVS*,2H**,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(5.,-9.,2.,5.,0.,-0.0)
C...  CALL QIKPLT(TP,P6,IP,8H*LOG(T)*,16H*CUM NO J EVALS*,2H**,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(5.,-9.,2.,5.,0.,1.0)
C...  CALL QIKPLT(TP,P7,IP,8H*LOG(T)*,16H*INC NO J EVALS*,2H**,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(5.,-9.,2.,5.,0.,1.0)
C...  CALL QIKPLT(TP,P8,IP,8H*LOG(T)*,15H*CURRENT ORDER*,2H**,1)
C...
C...  VAX PLOTTING AT SSCL
C...
C...  AT THE END OF THE CURRENT RUN, CREATE AN OUTPUT FILE FOR TOP DRAWER 
C...  PLOTTING
         OPEN(4,FILE='T.TOP',STATUS='NEW')
         NCASE=1
C...
C...  WRITE TOP DRAWER FILE FOR PLOTTING Y1 AND Y2
      WRITE(4,17)
17       FORMAT(' SET LIMITS X FROM -9 TO 1 Y FROM 0 TO 2',/,
     1          ' SET FONT DUPLEX')
      DO 13 IC=1,2
         WRITE(4,14)(TP(I),YP(IC,I),I=1,IP)
14       FORMAT(2F10.3)
         WRITE(4,16)
16       FORMAT(' JOIN 1')
13    CONTINUE
      WRITE(4,15)NCASE
15    FORMAT(
     1' Title 5.0 9.5 "Figure ',I2,'a: Y1, Y2 vs log(t)"'
     2,/,' TITLE LEFT "Y1, Y2"'
     3,/,' TITLE BOTTOM "log(t)"'
     4,/,' TITLE 4.5 0.75 '
     5   '"Y1, Y2 vs log(t), t = 1.0E-09,..., 1.0E+01"')
      WRITE(4,18)
18    FORMAT(' NEW FRAME')
C...
C...  WRITE TOP DRAWER FILE FOR PLOTTING LAST STEP SIZE
      WRITE(4,27)
27       FORMAT(' SET LIMITS X FROM -9 TO 1 Y FROM 0 TO 0.2',/,
     1          ' SET FONT DUPLEX')
      DO 23 IC=1,1
         WRITE(4,24)(TP(I),P1(I),I=1,IP)
24       FORMAT(2F10.3)
         WRITE(4,26)
26       FORMAT(' JOIN 1')
23    CONTINUE
      WRITE(4,25)NCASE
25    FORMAT(
     1' Title 4.0 9.5 "Figure ',I2,'b: HU (last step size) vs log(t)"'
     2,/,' TITLE LEFT "HU"'
     3,/,' TITLE BOTTOM "log(t)"'
     4,/,' TITLE 4.5 0.75 '
     5   '"HU vs log(t), t = 1.0E-09,..., 1.0E+01"')
      WRITE(4,18)
C...
C...  WRITE TOP DRAWER FILE FOR PLOTTING TOTAL STEPS
      WRITE(4,37)
37       FORMAT(' SET LIMITS X FROM -9 TO 1 Y FROM 0 TO 200',/,
     1          ' SET FONT DUPLEX')
      DO 33 IC=1,1
         WRITE(4,34)(TP(I),P2(I),I=1,IP)
34       FORMAT(2F10.3)
         WRITE(4,36)
36       FORMAT(' JOIN 1')
33    CONTINUE
      WRITE(4,35)NCASE
35    FORMAT(
     1' Title 4.0 9.5 "Figure ',I2,'c: NST (total steps) vs log(t)"'
     2,/,' TITLE LEFT "NST"'
     3,/,' TITLE BOTTOM "log(t)"'
     4,/,' TITLE 4.5 0.75 '
     5   '"NST vs log(t), t = 1.0E-09,..., 1.0E+01"')
      WRITE(4,18)
C...
C...  WRITE TOP DRAWER FILE FOR PLOTTING INCREMENTAL STEPS
      WRITE(4,47)
47       FORMAT(' SET LIMITS X FROM -9 TO 1 Y FROM 0 TO 50',/,
     1          ' SET FONT DUPLEX')
      DO 43 IC=1,1
         WRITE(4,44)(TP(I),P3(I),I=1,IP)
44       FORMAT(2F10.3)
         WRITE(4,46)
46       FORMAT(' JOIN 1')
43    CONTINUE
      WRITE(4,45)NCASE
45    FORMAT(
     1' Title 4.0 9.5 "Figure ',I2,'d: Incremental steps vs log(t)"'
     2,/,' TITLE LEFT "Inc NST"'
     3,/,' TITLE BOTTOM "log(t)"'
     4,/,' TITLE 4.5 0.75 '
     5   '"Inc NST vs log(t), t = 1.0E-09,..., 1.0E+01"')
      WRITE(4,18)
C...
C...  WRITE TOP DRAWER FILE FOR PLOTTING FUNCTION EVALUATIONS
      WRITE(4,57)
57       FORMAT(' SET LIMITS X FROM -9 TO 1 Y FROM 0 TO 300',/,
     1          ' SET FONT DUPLEX')
      DO 53 IC=1,1
         WRITE(4,54)(TP(I),P4(I),I=1,IP)
54       FORMAT(2F10.3)
         WRITE(4,56)
56       FORMAT(' JOIN 1')
53    CONTINUE
      WRITE(4,55)NCASE
55    FORMAT(
     1' Title 4.0 9.5 "Figure ',I2,'e: NFE (func eval) vs log(t)"'
     2,/,' TITLE LEFT "NFE"'
     3,/,' TITLE BOTTOM "log(t)"'
     4,/,' TITLE 4.5 0.75 '
     5   '"NFE vs log(t), t = 1.0E-09,..., 1.0E+01"')
      WRITE(4,18)
C...
C...  WRITE TOP DRAWER FILE FOR PLOTTING INCREMENTAL FUNCTION EVALUATIONS
      WRITE(4,67)
67       FORMAT(' SET LIMITS X FROM -9 TO 1 Y FROM 0 TO 50',/,
     1          ' SET FONT DUPLEX')
      DO 63 IC=1,1
         WRITE(4,64)(TP(I),P5(I),I=1,IP)
64       FORMAT(2F10.3)
         WRITE(4,66)
66       FORMAT(' JOIN 1')
63    CONTINUE
      WRITE(4,65)NCASE
65    FORMAT(
     1' Title 4.0 9.5 "Figure ',I2,'f: Inc function eval vs log(t)"'
     2,/,' TITLE LEFT "Inc NFE"'
     3,/,' TITLE BOTTOM "log(t)"'
     4,/,' TITLE 4.5 0.75 '
     5   '"Inc NFE vs log(t), t = 1.0E-09,..., 1.0E+01"')
      WRITE(4,18)
C...
C...  WRITE TOP DRAWER FILE FOR PLOTTING JACOBIAN EVALUATIONS
      WRITE(4,77)
77       FORMAT(' SET LIMITS X FROM -9 TO 1 Y FROM 0 TO 30',/,
     1          ' SET FONT DUPLEX')
      DO 73 IC=1,1
         WRITE(4,74)(TP(I),P6(I),I=1,IP)
74       FORMAT(2F10.3)
         WRITE(4,76)
76       FORMAT(' JOIN 1')
73    CONTINUE
      WRITE(4,75)NCASE
75    FORMAT(
     1' Title 4.0 9.5 "Figure ',I2,'g: NJE (Jacobian eval) vs log(t)"'
     2,/,' TITLE LEFT "NJE"'
     3,/,' TITLE BOTTOM "log(t)"'
     4,/,' TITLE 4.5 0.75 '
     5   '"NJE vs log(t), t = 1.0E-09,..., 1.0E+01"')
      WRITE(4,18)
C...
C...  WRITE TOP DRAWER FILE FOR PLOTTING INCREMENTAL JACOBIAN EVALUATIONS
      WRITE(4,87)
87       FORMAT(' SET LIMITS X FROM -9 TO 1 Y FROM 0 TO 3',/,
     1          ' SET FONT DUPLEX')
      DO 83 IC=1,1
         WRITE(4,84)(TP(I),P7(I),I=1,IP)
84       FORMAT(2F10.3)
         WRITE(4,86)
86       FORMAT(' JOIN 1')
83    CONTINUE
      WRITE(4,85)NCASE
85    FORMAT(
     1' Title 4.0 9.5 "Figure ',I2,'h: Inc Jacobian eval vs log(t)"'
     2,/,' TITLE LEFT "Inc NJE"'
     3,/,' TITLE BOTTOM "log(t)"'
     4,/,' TITLE 4.5 0.75 '
     5   '"Inc NJE vs log(t), t = 1.0E-09,..., 1.0E+01"')
      WRITE(4,18)
C...
C...  WRITE TOP DRAWER FILE FOR PLOTTING LAST ORDER
      WRITE(4,97)
97       FORMAT(' SET LIMITS X FROM -9 TO 1 Y FROM 0 TO 6',/,
     1          ' SET FONT DUPLEX')
      DO 93 IC=1,1
         WRITE(4,94)(TP(I),P8(I),I=1,IP)
94       FORMAT(2F10.3)
         WRITE(4,96)
96       FORMAT(' JOIN 1')
93    CONTINUE
      WRITE(4,95)NCASE
95    FORMAT(
     1' Title 4.0 9.5 "Figure ',I2,'i: NQU (last order) vs log(t)"'
     2,/,' TITLE LEFT "NQU"'
     3,/,' TITLE BOTTOM "log(t)"'
     4,/,' TITLE 4.5 0.75 '
     5   '"NQU vs log(t), t = 1.0E-09,..., 1.0E+01"')
C...  ******************************************************************
      RETURN
      END
