*APP MATH51
      PROGRAM IVPAGS
C...
C...  AN INTERFACE BETWEEN DSS/2 APPLICATIONS AND IMSL ROUTINE IVPAG
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 IMSL IVPAG (AN IMPLEMENTATION OF THE ADAMS' METHOD FOR NONSTIFF
C...  ODES, AND GEAR'S BDF METHOD FOR STIFF ODES)
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...  COMMON AREA TO LINK THIS MAIN PROGRAM WITH THE USER-SUPPLIED
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 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 IVPAG
      COMMON/S/ PARAM(50),   YV(250),    A(1,1)
C...
C...  EXTERNAL THE ROUTINES CALLED BY IVPAG
      EXTERNAL FCN, FCNJ
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 INPUT AND OUTPUT FILES
      OPEN(NI,FILE='DATA')
      OPEN(NO,FILE='OUTPUT')
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...  SET THE INITIAL CONDITIONS FOR SUBROUTINE IVPAG
      DO 5 I=1,NEQN
      YV(I)=Y(I)
5     CONTINUE
C...
C...  SET THE PARAMETERS FOR SUBROUTINE IVPAG
      TV=T0
      IDO=1
      TOL=1.0E-06
C...
C...  SET ARRAY PARAM TO THE DEFAULTS.  SUBROUTINE IVPAG HAS MANY
C...  OPTIONS AS SET BY THE ELEMENTS OF ARRAY PARAM WHICH CONTROL
C...  THE NUMERICAL INTEGRATION.  THE SELECTION OF THESE OPTIONS
C...  WILL, IN GENERAL, DETERMINE THE PERFORMANCE OF IVPAG, AND
C...  MAY, IN FACT, DETERMINE IF A SOLUTION CAN BE COMPUTED. THE
C...  TEST PROBLEM DEFINED IN SUBROUTINES INITAL, DERV AND PRINT
C...  IS MODEST IN SIZE AND COMPUTATIONAL REQUIREMENTS SO THAT MOST
C...  OF THE DEFAULTS FOR ARRAY PARAM CAN BE USED.  THE DETAILS OF THE
C...  OPTIONS ARE GIVEN IN THE IMSL DOCUMENTATION FOR IVPAG.
      DO 7 I=1,50
      PARAM(I)=0.
7     CONTINUE
C...
C...  SINCE THE PROBLEM IS STIFF, SELECT THE STIFF OPTION OF IVPAG WITH
C...  FINITE DIFFERENCE APPROXIMATION OF THE JACOBIAN MATRIX
      PARAM(12)=2.
      PARAM(13)=2.
C...
C...  CALL IVPAG TO COVER ONE PRINT INTERVAL
      TEND=TV+TP
4     CALL IVPAG(IDO,NEQN,FCN,FCNJ,A,TV,TEND,TOL,PARAM,YV)
C...
C...  TEST FOR AN ERROR CONDITION
      IF((IERROR.EQ.1).AND.(PARAM(8).NE.0.).AND.(IDO.EQ.6))THEN
C...
C...     PRINT A MESSAGE INDICATING AN ERROR CONDITION
         WRITE(NO,1004)PARAM(8),IDO
C...
C...     GO ON TO THE NEXT RUN
         GO TO 1
      END IF
C...
C...  ONE PRINT INTERVAL IS COMPLETE, SO PRINT THE SOLUTION
      T=TV
      DO 6 I=1,NEQN
      Y(I)=YV(I)
6     CONTINUE
      CALL PRINT(NI,NO)
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
      TEND=TV*(1.0D+01)**(1.0D+00/3.0D+00)

C...  CHECK FOR THE END OF THE RUN
      IF(TV.LT.(TF-0.5*TP))GO TO 4
C...
C...  THE CURRENT RUN IS COMPLETE, SO PRINT THE COMPUTATAIONAL
C...  STATISTICS AND GO ON TO THE NEXT RUN
      WRITE(NO,1005)(PARAM(I),I=32,36)
      GO TO 1
C...
C...  *****************************************************************
C...
C...  FORMATS
C...
1000  FORMAT(20A4)
1001  FORMAT(3E10.0)
1002  FORMAT(4I5,2X,A3,E10.0)
1003  FORMAT(1H1,
     1 ' RUN NO. - ',I3,2X,20A4,//,
     2 ' INITIAL T - ',E10.3,//,
     3 '   FINAL T - ',E10.3,//,
     4 '   PRINT T - ',E10.3,//,
     5 ' NUMBER OF DIFFERENTIAL EQUATIONS - ',I3,//,
     6 ' PRINT INTERVAL/MINIMUM INTEGRATION INTERVAL - ',I5,//,
     7 ' INTEGRATION ALGORITHM - ',I2,' - IVPAG ',//,
     8 ' INTEGRATION ERROR MESSAGES - ',I1,//,
     9 ' ERROR CRITERION - ',A3,//,
     A ' MAXIMUM INTEGRATION ERROR - ',E10.3,//,
     B 1H1)
1004  FORMAT(1H ,//,' PARAM(8) = ',E10.3,3X,' IDO = ',I2,//,
     1 ' INDICATING AN INTEGRATION ERROR, SO THE CURRENT RUN IS  ',/,
     2 ' TERMINATED.  PLEASE REFER TO THE IMSL DOCUMENTATION FOR ',/,
     3 ' SUBROUTINE IVPAG FOR AN EXPLANATION OF THESE ERROR INDI-',/,
     4 ' CATORS.                                                 ')
1005  FORMAT(//,' COMPUTATIONAL STATISTICS',/,
     1          ' MINIMUM STEP SIZE = ',E10.3,/,
     2          ' MAXIMUM STEP SIZE = ',E10.3,/,
     3          ' TOTAL NUMBER OF STEPS = ',E10.3,/,
     4          ' NUMBER OF FUNCTION EVALUATIONS = ',E10.3,/,
     5          ' NUMBER OF JACOBIAN EVALUATIONS = ',E10.3)
      END
      SUBROUTINE FCN(NEQN,TV,YV,YDOT)
C...
C...  SUBROUTINE FCN IS AN INTERFACE ROUTINE BETWEEN SUBROUTINES IVPAG
C...  AND DERV
C...
C...  DSS/2 COMMON AREA
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(1)
     2      /F/      F(1)
C...
C...  VARIABLE DIMENSION THE DEPENDENT, DERIVATIVE VECTORS
      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 IVPAG
      DO 2 I=1,NEQN
      YDOT(I)=F(I)
2     CONTINUE
      RETURN
      END
      SUBROUTINE FCNJ(NEQN,T,Y,DYPDY)
C...
C...  SUBROUTINE FCNJ COMPUTES THE JACOBIAN MATRIX OF THE ODE SYSTEM.
C...  FOR THIS APPLICATION, FCNJ IS NOT CALLED BECAUSE
C...
C...     (1)  THE NONSTIFF OPTION, ADAMS' METHOD, IS SELECTED, I.E.,
C...
C...          METH = 1, PARAM(12) = 1.
C...
C...     (2)  THE STIFF OPTION, GEAR'S (BDF) METHOD, IS SELECTED, WITH
C...          CHORD METHOD AND INTERNALLY GENERATED FINITE-DIFFERENCE
C...          JACOBIAN, I.E.,
C...
C...          METH = 2, MITER = 2, PARAM(12) = 2., PARAM(13) = 2.
C...
      REAL Y(NEQN), DYPDY(*)
      RETURN
      END
      SUBROUTINE INITAL
C...
C...  TWO ODES WITH A STIFFNESS RATIO OF 1.0E+06 INTEGRATED BY IMSL
C...  ROUTINE IVPAG USING DSS/2 FORMAT - PRINTING OF THE NUMERICAL
C...  AND EXACT SOLUTIONS
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
      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.E+00
      Y2=2.E+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.5E+00
      B=499999.5E+00
      RETURN
      END
      SUBROUTINE DERV
      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)
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/    Y1,    Y2
     2      /F/ DY1DT, DY2DT
     3      /C/     A,     B
C...
C...  COMMON/S/ ACCESSES THE COMPUTATIONAL STATISTICS OF IVPAG
      COMMON/S/ PARAM(50)
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION    TP(31),
     1             P1(31),    P2(31),    P3(31),    P4(31),    P5(31)
C...
C...  INITIALIZE A COUNTER FOR THE PLOTTED STATISTICS
      DATA IP/0/
C...
C...  PRINT A HEADING FOR THE COMPUTATIONAL STATISTICS
      IF(IP.EQ.0)WRITE(NO,2)
2     FORMAT(1H ,//,9X,'T',7X,'MIN H',7X,'MAX H',7X,'STEPS',
     1       6X,'F EVAL',6X,'J EVAL')
C...
C...  EVALUATE THE EXACT (ANALYTICAL) SOLUTION, EQUATIONS (11) AND (12)
C...  WITH L1 = -(A + B), L2 = -(A - B), FOR COMPARISON WITH THE NUMERI-
C...  CAL SOLUTION
      AMB=A-B
      IF(AMB.EQ.0.E+00)AMB=1.E+00
      EXP1=EXP(-AMB*T)
      EXP2=0.E+00
      IF(((A+B)*T).LT.50.E+00)EXP2=EXP(-(A+B)*T)
      Y1E=EXP1-EXP2
      Y2E=EXP1+EXP2
C...
C...  COMPUTE THE DIFFERENCES BETWEEN THE NUMERICAL AND EXACT SOLUTIONS
      Y1D=Y1-Y1E
      Y2D=Y2-Y2E
C...
C...  INITIALIZE THREE VARIABLES FOR COMPUTING THE NET NUMBER OF STEPS,
C...  FUNCTION EVALUATIONS AND JACOBIAN MATRIX EVALUATIONS
      WRITE(NO,1)T,(PARAM(I),I=32,36)
1     FORMAT(E10.3,5E12.4)
C...
C...  STORE THE COMPUTATIONAL STATISTICS FOR SUBSEQUENT PLOTTING
      IP=IP+1
      TP(IP)=ALOG10(T)
      P1(IP)=PARAM(32)
      P2(IP)=PARAM(33)
      P3(IP)=PARAM(34)
      P4(IP)=PARAM(35)
      P5(IP)=PARAM(36)
C...
C...  PLOT THE 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...  CALL QIKSAX(3,3)
C...  CALL QIKSET(5.,-9.,2.,5.,0.,-0.0)
C...  CALL QIKPLT(TP,P1,IP,8H*LOG(T)*,7H*MIN H*,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)*,7H*MAX H*,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)*,10H*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)*,8H*F EVAL*,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)*,8H*J EVAL*,2H**,1)
C...
C...  ******************************************************************
      RETURN
      END
LAPIDUS, ET AL, NUM METH DIFF SYS, P107
    0.D+00    1.D+01    1.D-09
    2 1000    2    1  REL 0.00001
END OF RUNS
*APP MATH52
      PROGRAM IVPAGS
C...
C...  AN INTERFACE BETWEEN DSS/2 APPLICATIONS AND IMSL ROUTINE IVPAG
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 IMSL IVPAG (AN IMPLEMENTATION OF THE ADAMS' METHOD FOR NONSTIFF
C...  ODES, AND GEAR'S BDF METHOD FOR STIFF ODES)
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...  COMMON AREA TO LINK THIS MAIN PROGRAM WITH THE USER-SUPPLIED
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 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 IVPAG
      PARAMETER(MXPARM=50,  NEQ=250)
      REAL  PARAM(MXPARM),   YV(NEQ),    A(1,1)
C...
C...  EXTERNAL THE ROUTINES CALLED BY IVPAG
      EXTERNAL FCN, FCNJ
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 INPUT AND OUTPUT FILES
      OPEN(NI,FILE='DATA')
      OPEN(NO,FILE='OUTPUT')
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 IVPAG
      DO 5 I=1,NEQN
      YV(I)=Y(I)
5     CONTINUE
C...
C...  SET THE PARAMETERS FOR SUBROUTINE IVPAG
      TV=T0
      IDO=1
      TOL=1.0E-06
C...
C...  SET ARRAY PARAM TO THE DEFAULTS.  SUBROUTINE IVPAG HAS MANY
C...  OPTIONS AS SET BY THE ELEMENTS OF ARRAY PARAM WHICH CONTROL
C...  THE NUMERICAL INTEGRATION.  THE SELECTION OF THESE OPTIONS
C...  WILL, IN GENERAL, DETERMINE THE PERFORMANCE OF IVPAG, AND
C...  MAY, IN FACT, DETERMINE IF A SOLUTION CAN BE COMPUTED. THE
C...  TEST PROBLEM DEFINED IN SUBROUTINES INITAL, DERV AND PRINT
C...  IS MODEST IN SIZE AND COMPUTATIONAL REQUIREMENTS SO THAT MOST
C...  OF THE DEFAULTS FOR ARRAY PARAM CAN BE USED.  THE DETAILS OF THE
C...  OPTIONS ARE GIVEN IN THE IMSL DOCUMENTATION FOR IVPAG.
      DO 7 I=1,50
      PARAM(I)=0.
7     CONTINUE
C...
C...  SINCE THE PROBLEM IS STIFF, SELECT THE STIFF OPTION OF IVPAG WITH
C...  FINITE DIFFERENCE APPROXIMATION OF THE JACOBIAN MATRIX
      PARAM(12)=2.
      PARAM(13)=2.
C...
C...  CALL IVPAG TO COVER ONE PRINT INTERVAL
4     TEND=TV*10.
      CALL IVPAG(IDO,NEQN,FCN,FCNJ,A,TV,TEND,TOL,PARAM,YV)
C...
C...  TEST FOR AN ERROR CONDITION
      IF((IERROR.EQ.1).AND.(PARAM(8).NE.0.).AND.(IDO.EQ.6))THEN
C...
C...     PRINT A MESSAGE INDICATING AN ERROR CONDITION
         WRITE(NO,1004)PARAM(8),IDO
C...
C...     GO ON TO THE NEXT RUN
         GO TO 1
      END IF
C...
C...  ONE PRINT INTERVAL IS COMPLETE, SO PRINT THE SOLUTION
      T=TV
      DO 6 I=1,NEQN
      Y(I)=YV(I)
6     CONTINUE
      CALL PRINT(NI,NO)
C...
C...  CHECK FOR A RUN TERMINATION
      IF(NSTOP.NE.0)GO TO 1
C...
C...  CHECK FOR THE END OF THE RUN
      IF(TV.LT.(TF-0.5*TP))GO TO 4
C...
C...  THE CURRENT RUN IS COMPLETE, SO PRINT THE COMPUTATAIONAL
C...  STATISTICS AND GO ON TO THE NEXT RUN
      WRITE(NO,1005)(PARAM(I),I=32,36)
      GO TO 1
C...
C...  *****************************************************************
C...
C...  FORMATS
C...
1000  FORMAT(20A4)
1001  FORMAT(3E10.0)
1002  FORMAT(4I5,2X,A3,E10.0)
1003  FORMAT(1H1,
     1 ' RUN NO. - ',I3,2X,20A4,//,
     2 ' INITIAL T - ',E10.3,//,
     3 '   FINAL T - ',E10.3,//,
     4 '   PRINT T - ',E10.3,//,
     5 ' NUMBER OF DIFFERENTIAL EQUATIONS - ',I3,//,
     6 ' PRINT INTERVAL/MINIMUM INTEGRATION INTERVAL - ',I5,//,
     7 ' INTEGRATION ALGORITHM - ',I2,' - IVPAG ',//,
     8 ' INTEGRATION ERROR MESSAGES - ',I1,//,
     9 ' ERROR CRITERION - ',A3,//,
     A ' MAXIMUM INTEGRATION ERROR - ',E10.3,//,
     B 1H1)
1004  FORMAT(1H ,//,' PARAM(8) = ',E10.3,3X,' IDO = ',I2,//,
     1 ' INDICATING AN INTEGRATION ERROR, SO THE CURRENT RUN IS  ',/,
     2 ' TERMINATED.  PLEASE REFER TO THE IMSL DOCUMENTATION FOR ',/,
     3 ' SUBROUTINE IVPAG FOR AN EXPLANATION OF THESE ERROR INDI-',/,
     4 ' CATORS.                                                 ')
1005  FORMAT(//,' COMPUTATIONAL STATISTICS',/,
     1          ' MINIMUM STEP SIZE = ',E10.3,/,
     2          ' MAXIMUM STEP SIZE = ',E10.3,/,
     3          ' TOTAL NUMBER OF STEPS = ',E10.3,/,
     4          ' NUMBER OF FUNCTION EVALUATIONS = ',E10.3,/,
     5          ' NUMBER OF JACOBIAN EVALUATIONS = ',E10.3)
      END
      SUBROUTINE FCN(NEQN,TV,YV,YDOT)
C...
C...  SUBROUTINE FCN IS AN INTERFACE ROUTINE BETWEEN SUBROUTINES IVPAG
C...  AND DERV
C...
C...  DSS/2 COMMON AREA
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(1)
     2      /F/      F(1)
C...
C...  VARIABLE DIMENSION THE DEPENDENT, DERIVATIVE VECTORS
      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 IVPAG
      DO 2 I=1,NEQN
      YDOT(I)=F(I)
2     CONTINUE
      RETURN
      END
      SUBROUTINE FCNJ(NEQN,T,Y,DYPDY)
C...
C...  SUBROUTINE FCNJ COMPUTES THE JACOBIAN MATRIX OF THE ODE SYSTEM.
C...  FOR THIS APPLICATION, FCNJ IS NOT CALLED BECAUSE
C...
C...     (1)  THE NONSTIFF OPTION, ADAMS' METHOD, IS SELECTED, I.E.,
C...
C...          METH = 1, PARAM(12) = 1.
C...
C...     (2)  THE STIFF OPTION, GEAR'S (BDF) METHOD, IS SELECTED, WITH
C...          CHORD METHOD AND INTERNALLY GENERATED FINITE-DIFFERENCE
C...          JACOBIAN, I.E.,
C...
C...          METH = 2, MITER = 2, PARAM(12) = 2., PARAM(13) = 2.
C...
      REAL Y(NEQN), DYPDY(*)
      RETURN
      END
      SUBROUTINE INITAL
C...
C...  TWO ODES WITH A STIFFNESS RATIO OF 1.0E+06 INTEGRATED BY IMSL
C...  ROUTINE IVPAG USING DSS/2 FORMAT - PRINTING OF THE COMPUTATIONAL
C...  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...
      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.E+00
      Y2=2.E+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.5E+00
      B=499999.5E+00
      RETURN
      END
      SUBROUTINE DERV
      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)
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/    Y1,    Y2
     2      /F/ DY1DT, DY2DT
     3      /C/     A,     B
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IF(T.EQ.0.E+00)WRITE(NO,2)
2     FORMAT(1H ,//,9X,1HT,
     1 13X,2HY1,12X,3HY1E,12X,3HY1D,/,
     2 23X,2HY2,12X,3HY2E,12X,3HY2D)
C...
C...  EVALUATE THE EXACT (ANALYTICAL) SOLUTION, EQUATIONS (11) AND (12)
C...  WITH L1 = -(A + B), L2 = -(A - B), FOR COMPARISON WITH THE NUMERI-
C...  CAL SOLUTION
      AMB=A-B
      IF(AMB.EQ.0.E+00)AMB=1.E+00
      EXP1=EXP(-AMB*T)
      EXP2=0.E+00
      IF(((A+B)*T).LT.50.E+00)EXP2=EXP(-(A+B)*T)
      Y1E=EXP1-EXP2
      Y2E=EXP1+EXP2
C...
C...  COMPUTE THE DIFFERENCES BETWEEN THE NUMERICAL AND EXACT SOLUTIONS
      Y1D=Y1-Y1E
      Y2D=Y2-Y2E
C...
C...  PRINT THE INDEPENDENT VARIABLE, NUMERICAL AND EXACT SOLUTIONS, AND
C...  THE DIFFERENCES BETWEEN THE TWO SOLUTIONS
      WRITE(NO,1)T,Y1,Y1E,Y1D,Y2,Y2E,Y2D
1     FORMAT(1H ,/,E10.3,3E15.4,/,10X,3E15.4)
      RETURN
      END
LAPIDUS, ET AL, NUM METH DIFF SYS, P107
    1.E-10    1.E+01    1.
    2 1000    2    1  REL 0.00001
END OF RUNS
*APP MATH53
       SUBROUTINE INITAL
 C...
 C...  SOLUTION OF TWO FOURTH-ORDER PDES
 C...
 C...  CONSIDER THE TWO PROBLEMS
 C...
 C...  (1)
 C...
 C...     U  = -U
 C...      T     XXXX
 C...
 C...     U(X,0) = COS(X/A)
 C...
 C...     UX(-PI*A,T) = UX(PI*A,T) = 0
 C...
 C...     UXXX(-PI*A,T) = UXXX(PI*A,T) = 0
 C...
 C...  WHICH HAS THE EXACT SOLUTION
 C...
 C...     U(X,T) = COS(X/A)*EXP((-1/A**4)*T)
 C...
 C...  (2)
 C...
 C...     U   = -U      (FOURTH ORDER BEAM EQUATION)
 C...      TT     XXXX
 C...
 C...     U(X,0) = COS(X/A), UT(X,0) = 0
 C...
 C...     UX(-PI*A,T) = UX(PI*A,T) = 0
 C...
 C...     UXXX(-PI*A,T) = UXXX(PI*A,T) = 0
 C...
 C...  WHICH HAS THE EXACT SOLUTION
 C...
 C...     U(X,T) = COS(X/A)*COS((1/A**2)*T)
 C...
 C...  NOTE - THE NUMERICAL SOLUTION FOR PROBLEM (2) WHICH FOLLOWS NEEDS
 C...  SOME MORE WORK WITH REGARD TO ACCURACY AND STABILITY, E.G., A STIFF
 C...  INTEGRATOR TO IMPROVE THE STABILITY OF THE CALCULATIONS, THEREBY
 C...  PERMITTING A LARGER TIME STEP (AND FINAL VALUE OF T).
 C...
       COMMON/T/            T,        NSTOP,        NORUN
 C...
 C...  INITIAL CONDITIONS FOR PROBLEM 1
       IF(NORUN.EQ.1)CALL INIT1
 C...
 C...  INITIAL CONDITIONS FOR PROBLEM 2
       IF(NORUN.EQ.2)CALL INIT2
       RETURN
       END
       SUBROUTINE DERV
       COMMON/T/            T,        NSTOP,        NORUN
 C...
 C...  DERIVATIVES FOR PROBLEM 1
       IF(NORUN.EQ.1)CALL DERV1
 C...
 C...  DERIVATIVES FOR PROBLEM 2
       IF(NORUN.EQ.2)CALL DERV2
       RETURN
       END
1      SUBROUTINE PRINT(NI,NO)
       COMMON/T/            T,        NSTOP,        NORUN
 C...
 C...  PRINTING FOR PROBLEM 1
       IF(NORUN.EQ.1)CALL PRINT1(NI,NO)
 C...
 C...  PRINTING FOR PROBLEM 2
       IF(NORUN.EQ.2)CALL PRINT2(NI,NO)
       RETURN
       END
       SUBROUTINE INIT1
       COMMON/T/            T,        NSTOP,        NORUN
      1      /Y/    U(-10:10)
      2      /F/   UT(-10:10)
      3      /S/   UX(-10:10),  UXX(-10:10), UXXX(-10:10),UXXXX(-10:10)
      4      /C/            A,           PI,    X(-10:10),            N,
      5                    IP
 C...
 C...  MODEL PARAMETERS
       A=1.0
       PI=ACOS(-1.0)
 C...
 C...  SET THE INITIAL CONDITIONS
       N=10
       DO 1 I=-N,N
       X(I)=A*PI*FLOAT(I)/FLOAT(N)
       U(I)=COS(X(I)/A)
 1     CONTINUE
 C...
 C...  INITIAL DERIVATIVES
       CALL DERV1
       IP=0
       RETURN
       END
       SUBROUTINE DERV1
       COMMON/T/            T,        NSTOP,        NORUN
      1      /Y/    U(-10:10)
      2      /F/   UT(-10:10)
      3      /S/   UX(-10:10),  UXX(-10:10), UXXX(-10:10),UXXXX(-10:10)
      4      /C/            A,           PI,    X(-10:10),            N,
      5                    IP
 C...
 C...  FIRST DERIVATIVE
       XL=-A*PI
       XU= A*PI
       NT=2*N+1
       CALL DSS006(XL,XU,NT,U,UX)
 C...
 C...  FIRST AND SECOND BOUNDARY CONDITIONS
       UX( N)=0.
       UX(-N)=0.
 C...
 C...  SECOND DERIVATIVE
       CALL DSS006(XL,XU,NT,UX,UXX)
 C...
 C...  THIRD DERIVATIVE
       CALL DSS006(XL,XU,NT,UXX,UXXX)
 C...
 C...  THIRD AND FOURTH BOUNDARY CONDITIONS
       UXXX( N)=0.
1      UXXX(-N)=0.
 C...
 C...  FOURTH DERIVATIVE
       CALL DSS006(XL,XU,NT,UXXX,UXXXX)
 C...
 C...  PDE
       DO 1 I=-N,N
       UT(I)=-UXXXX(I)
 1     CONTINUE
       RETURN
       END
       SUBROUTINE PRINT1(NI,NO)
       COMMON/T/            T,        NSTOP,        NORUN
      1      /Y/    U(-10:10)
      2      /F/   UT(-10:10)
      3      /S/   UX(-10:10),  UXX(-10:10), UXXX(-10:10),UXXXX(-10:10)
      4      /C/            A,           PI,    X(-10:10),            N,
      5                    IP
 C...
 C...  DIMENSION THE ARRAYS FOR THE ERRORS
       REAL         E(-10:10),   ET(-10:10),
      1            EX(-10:10),  EXX(-10:10), EXXX(-10:10),EXXXX(-10:10)
 C...
 C...  DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION
       REAL XP(21), UP(6,21)
 C...
 C...  PRINT THE NUMERICAL SOLUTION
       WRITE(NO,1)T,(X(I),I=-N,N),
      1                (U(I),I=-N,N),
      1               (UX(I),I=-N,N),
      1              (UXX(I),I=-N,N),
      1             (UXXX(I),I=-N,N),
      1            (UXXXX(I),I=-N,N),
      1               (UT(I),I=-N,N)
 1     FORMAT(//,' T = ',F6.1,/,
      1'     X',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'     U',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'    UX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'   UXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'  UXXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1' UXXXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'    UT',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//)
 C...
 C...  CALCULATE THE ERRORS IN THE NUMERICAL SOLUTION
       DO 2 I=-N,N
           SINE=SIN(X(I)/A)
           COSE=COS(X(I)/A)
          EXPON=EXP((-1.0/A**4)*T)
           E(I)=    U(I)-( 1.0     )*(COSE*EXPON)
          EX(I)=   UX(I)-(-1.0/A   )*(SINE*EXPON)
         EXX(I)=  UXX(I)-(-1.0/A**2)*(COSE*EXPON)
        EXXX(I)= UXXX(I)-(+1.0/A**3)*(SINE*EXPON)
       EXXXX(I)=UXXXX(I)-(+1.0/A**4)*(COSE*EXPON)
          ET(I)=   UT(I)-(-1.0/A**4)*(COSE*EXPON)
 2     CONTINUE
 C...
 C...  PRINT THE ERRORS
       WRITE(NO,3)
      1                (E(I),I=-N,N),
      1               (EX(I),I=-N,N),
1     1              (EXX(I),I=-N,N),
      1             (EXXX(I),I=-N,N),
      1            (EXXXX(I),I=-N,N),
      1               (ET(I),I=-N,N)
 3     FORMAT(//,
      1'     E',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'    EX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'   EXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'  EXXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1' EXXXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'    ET',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//)
 C...
 C...  STORE THE SOLUTION FOR PLOTTING
       IP=IP+1
       DO 4 I=-N,N
       XP(I+N+1)=X(I)
       UP(IP,I+N+1)=U(I)
 4     CONTINUE
 C...
 C...  PLOT THE SOLUTION
       IF(IP.LT.6)RETURN
       NT=2*N+1
       CALL TPLOTS(IP,NT,XP,UP)
       WRITE(NO,5)
 5     FORMAT(/,' U(X,T) VS X, PROBLEM 1, T = 0, 1, 2,..., 5')
       RETURN
       END
       SUBROUTINE INIT2
       COMMON/T/            T,        NSTOP,        NORUN
      1      /Y/    U(-10:10),   U1(-10:10)
      2      /F/   UT(-10:10),  U1T(-10:10)
      3      /S/   UX(-10:10),  UXX(-10:10), UXXX(-10:10),UXXXX(-10:10)
      4      /C/            A,           PI,    X(-10:10),            N,
      5                    IP
 C...
 C...  MODEL PARAMETERS
       A=1.0
       PI=ACOS(-1.0)
 C...
 C...  SET THE INITIAL CONDITIONS
       N=10
       DO 1 I=-N,N
       X(I)=A*PI*FLOAT(I)/FLOAT(N)
        U(I)=COS(X(I)/A)
       U1(I)=0.
 1     CONTINUE
 C...
 C...  INITIAL DERIVATIVES
       CALL DERV2
       IP=0
       RETURN
       END
       SUBROUTINE DERV2
       COMMON/T/            T,        NSTOP,        NORUN
      1      /Y/    U(-10:10),   U1(-10:10)
      2      /F/   UT(-10:10),  U1T(-10:10)
      3      /S/   UX(-10:10),  UXX(-10:10), UXXX(-10:10),UXXXX(-10:10)
      4      /C/            A,           PI,    X(-10:10),            N,
      5                    IP
 C...
1C...  FIRST DERIVATIVE
       XL=-A*PI
       XU= A*PI
       NT=2*N+1
       CALL DSS006(XL,XU,NT,U,UX)
 C...
 C...  FIRST AND SECOND BOUNDARY CONDITIONS
       UX( N)=0.
       UX(-N)=0.
 C...
 C...  SECOND DERIVATIVE
       CALL DSS006(XL,XU,NT,UX,UXX)
 C...
 C...  THIRD DERIVATIVE
       CALL DSS006(XL,XU,NT,UXX,UXXX)
 C...
 C...  THIRD AND FOURTH BOUNDARY CONDITIONS
       UXXX( N)=0.
       UXXX(-N)=0.
 C...
 C...  FOURTH DERIVATIVE
       CALL DSS006(XL,XU,NT,UXXX,UXXXX)
 C...
 C...  PDE
       DO 1 I=-N,N
        UT(I)= U1(I)
       U1T(I)=-UXXXX(I)
 1     CONTINUE
       RETURN
       END
       SUBROUTINE PRINT2(NI,NO)
       COMMON/T/            T,        NSTOP,        NORUN
      1      /Y/    U(-10:10),   U1(-10:10)
      2      /F/   UT(-10:10),  U1T(-10:10)
      3      /S/   UX(-10:10),  UXX(-10:10), UXXX(-10:10),UXXXX(-10:10)
      4      /C/            A,           PI,    X(-10:10),            N,
      5                    IP
 C...
 C...  DIMENSION THE ARRAYS FOR THE ERRORS
       REAL         E(-10:10),   ET(-10:10),  E1T(-10:10),
      1            EX(-10:10),  EXX(-10:10), EXXX(-10:10),EXXXX(-10:10)
 C...
 C...  DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION
       REAL XP(21), UP(6,21)
 C...
 C...  PRINT THE NUMERICAL SOLUTION
       WRITE(NO,1)T,(X(I),I=-N,N),
      1                (U(I),I=-N,N),
      1               (UX(I),I=-N,N),
      1              (UXX(I),I=-N,N),
      1             (UXXX(I),I=-N,N),
      1            (UXXXX(I),I=-N,N),
      1               (UT(I),I=-N,N),
      1              (U1T(I),I=-N,N)
 1     FORMAT(//,' T = ',F6.1,/,
      1'     X',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'     U',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'    UX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'   UXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'  UXXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
1     1' UXXXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'    UT',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'   U1T',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//)
 C...
 C...  CALCULATE THE ERRORS IN THE NUMERICAL SOLUTION
       DO 2 I=-N,N
           SINX=SIN(X(I)/A)
           COSX=COS(X(I)/A)
           SINT=SIN(T/A**2)
           COST=COS(T/A**2)
           E(I)=    U(I)-( 1.0     )*(COSX*COST)
          EX(I)=   UX(I)-(-1.0/A   )*(SINX*COST)
         EXX(I)=  UXX(I)-(-1.0/A**2)*(COSX*COST)
        EXXX(I)= UXXX(I)-(+1.0/A**3)*(SINX*COST)
       EXXXX(I)=UXXXX(I)-(+1.0/A**4)*(COSX*COST)
          ET(I)=   UT(I)-(-1.0/A**2)*(COSX*SINT)
         E1T(I)=  U1T(I)-(-1.0/A**4)*(COSX*COST)
 2     CONTINUE
 C...
 C...  PRINT THE ERRORS
       WRITE(NO,3)
      1                (E(I),I=-N,N),
      1               (EX(I),I=-N,N),
      1              (EXX(I),I=-N,N),
      1             (EXXX(I),I=-N,N),
      1            (EXXXX(I),I=-N,N),
      1               (ET(I),I=-N,N),
      1              (E1T(I),I=-N,N)
 3     FORMAT(//,
      1'     E',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'    EX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'   EXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'  EXXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1' EXXXX',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'    ET',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//,
      1'   E1T',5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,5E12.3,/,6X,E12.3,//)
 C...
 C...  STORE THE SOLUTION FOR PLOTTING
       IP=IP+1
       DO 4 I=-N,N
       XP(I+N+1)=X(I)
       UP(IP,I+N+1)=U(I)
 4     CONTINUE
 C...
 C...  PLOT THE SOLUTION
       IF(IP.LT.6)RETURN
       NT=2*N+1
       CALL TPLOTS(IP,NT,XP,UP)
       WRITE(NO,5)
 5     FORMAT(/,' U(X,T) VS X, PROBLEM 2, T = 0, 1, 2,..., 5')
       RETURN
       END
FOURTH-ORDER PDES - PROBLEM 1
0.        5.0       1.0
   2199999    5    1  REL 0.001
FOURTH-ORDER PDES - PROBLEM 2
0.        1.0       0.2
   4299999    5    1  REL 0.001
END OF RUNS
*APP MATH54
      PROGRAM LSODE2
C...
C...  ADAPTIVE GRID SOLUTION OF THE ONE-DIMENSIONAL BURGERS EQUATION
C...  WITH BDF TIME INTEGRATION VIA LSODE
C...
C...  THE FOLLOWING PROGRAM CALLS THE USUAL USER-SUPPLIED SUBROUTINES
C...  INITAL, DERV AND PRINT, PLUS DATA, FOR THE NUMERICAL METHOD
C...  OF LINES SOLUTION OF THE ONE-DIMENSIONAL BURGERS EQUATION BY
C...  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...  ALSO, SUBROUTINE PRINT CALLS THE ADAPTIVE GRID SUBROUTINE ANUGB1
C...  TO REDFINE THE ADAPTIVE GRID AT EACH PRINT TIME.  ANUGB1 IN TURN
C...  CALLS A SERIES OF SUBORDINATE ROUTINES PROVIDED WITH DSS/2.
C...
C...  NOTE THAT THE FOLLOWING CODING IS IN DOUBLE PRECISION FORMAT
C...  AND THEREFORE REQUIRES THE DOUBLE PRECISION LSODE, ANUGB1 AND
C...  SUBORDINATE ROUTINES.
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 AREA TO PROVIDE THE INPUT/OUTPUT UNIT NUMBERS TO OTHER
C...  SUBROUTINES
      COMMON/IO/       NI,        NO
C...
C...  COMMON/SYSTM1/ IS ACCESSED SO THE NUMBER OF ODES DEFINED BY THE
C...  ADAPTIVE GRID, NEQN, CAN BE CHANGED AS THE GRID EVOLVES.  THE
C...  ADAPTIVE GRID IS REDEFINED VIA SUBROUTINE ANUGB1 CALLED FROM
C...  SUBROUTINE PRINT.  THEREFORE, THE ADAPTIVE GRID, AND THUS THE
C...  NUMBER OF ODES, NEQN, IS CHANGED AT THE END OF EACH PRINT IN-
C...  TERVAL.  THE NEW NUMBER OF ODES IS USED BY LSODE AS AN INPUT
C...  (SEE THE ADDITIONAL COMMENTS JUST BEFORE THE CALL TO LSODE BELOW)
      COMMON/SYSTM1/ T0,TF,TP,H,ERROR,NEQN
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 INPUT AND OUTPUT FILES
      OPEN(NI,FILE='DATA')
      OPEN(NO,FILE='OUTPUT')
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=1.0D-05
      ATOL=1.0D-03
      LRW=6600
      LIW=275
      IOPT=0
      ITASK=1
      ISTATE=1
      METH=2
      MITER=2
C...
C...  THE METHOD FLAG, MF, IS COMPUTED FROM METH AND MITER
      MF=10*METH+MITER
C...
C...  CALL SUBROUTINE LSODE TO COVER ONE PRINT INTERVAL.  NOTE THAT LSODE
C...  IS INITALIZED EACH TIME IT IS CALLED (ISTATE = 1) SO THAT THE ORDER
C...  OF THE INTEGRATOR IS THEREFORE RESET TO 1 (THE IMPLICIT EULER
C...  METHOD).  THIS OCCURS AT THE END OF EACH PRINT INTERVAL BECAUSE
C...  THE ADAPTIVE GRID IS REDEFINED BY A CALL TO SUBROUTINE ANUGB1 IN
C...  SUBROUTINE PRINT.  THIS REDEFINTION OF THE ADAPTIVE GRID, IN
C...  EFFECT, ALSO DEFINES A NEW ODE PROBLEM (THE NUMBER OF ODES, NEQN,
C...  CAN CHANGE THROUGH COMMON/SYSTM1/), AND THEREFORE LSODE MUST BE
C...  INITIALIZED FOR THE NEW ODE PROBLEM
4     ISTATE=1
      TOUT=TV+TP
      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...  SET THE INITIAL CONDITIONS FOR THE NEW SET OF ODES DEFINED BY THE
C...  CALL TO THE ADAPTIVE GRID ROUTINE ANUGB1 IN SUBROUTINE PRINT
      DO 7 I=1,NEQN
      YV(I)=Y(I)
7     CONTINUE
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...  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(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...  SOLUTION OF BURGERS' EQUATION BY SPLINE APPROXIMATION ON AN
C...  ADAPTIVE GRID
C...
C...  DOUBLE PRECISION
C...
C...  SUBROUTINE INITAL
C...
C...     (1)  DEFINES THE INITIAL ADAPTIVE GRID AS THE BASIC UNIFORM
C...          GRID VIA A CALL TO SUBROUTINE GRIDE
C...
C...     (2)  DEFINES THE INITIAL CONDITION ON THE ADAPTIVE GRID OF
C...          (1) BY A CALL TO FUNCTION PHI WHICH IMPLEMENTS THE EXACT
C...          SOLUTION TO BURGERS EQUATION
C...
C...  TYPE REAL VARIABLES AS DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  COMMON/SYSTM1/ IS ACCESSED SO THE NUMBER OF ODES DEFINED BY THE
C...  ADAPTIVE GRID, N, CAN BE CHANGED AS THE GRID EVOLVES
      COMMON/SYSTM1/ T0,TF,TP,H,ERROR,N
C...
C...  THE FOLLOWING COMMON AREAS ARE BRIEFLY EXPLAINED BELOW
C...
C...     /T/ /Y/ /F/    STANDARD DSS/2 COMMON AREAS DIMENSIONED FOR UP
C...                    TO 150 ODES
C...
C...     /B/            THE SOLUTION, CB, AND ITS SECOND DERIVATIVE,
C...                    CBXX, DEFINED ON THE BASIC GRID
C...
C...     /SV/           INDEPENDENT VARIABLE DEFINED ON THE ADAPTIVE
C...                    GRID
C...
C...     /SC/           CUBIC SPLINE COEFFICIENTS (LINEAR, QUADRATIC
C...                    AND CUBIC COEFFICIENTS RESPECTIVELY)
C...
C...     /C/            CONSTANTS PERTAINING TO THE PROBLEM
C...
C...     /AB1/          PARAMETERS WHICH CONTROL THE ADAPTIVE GRID
C...
      COMMON  /T/         T
     1        /Y/   CA(150)
     2        /F/   CT(150)
     3        /B/   CB(120), CBXX(120)
     4       /SV/   XA(150)
     5       /SC/   CX(150),  CXX(150), CXXX(150)
     6        /C/        XL,        XU,         V,       VIS,
     7                   SN,        NB,        IP,        JP
      COMMON/AB1/        DL,       XFL,       XFR,        NF,
     1                  DXF,   XB(120),  LXB(120),  NOB(120)
C...
C...  SET THE VELOCITY AND VISCOSITY OF BURGERS EQUATION
      V=1.0D+00
      VIS=0.003D+00
C...
C...  INITIALIZE THE COUNTERS FOR PRINTED AND PLOTTED OUTPUT USED
C...  IN SUBROUTINE PRINT, THE RUNNING SUM OF THE NUMBER OF ODES
C...  USED IN THE ADAPTIVE GRID
      IP=0
      JP=0
      SN=0.0D+00
C...
C...  DEFINE THE LENGTH OF THE X AXIS
      XL=0.0D+00
      XU=1.0D+00
C...
C...  DEFINE THE NUMBER OF BASIC GRID POINTS, NUMBER OF SUBINTERVALS
C...  BETWEEN THE BASIC GRID POINTS WHEN ADAPTIVE POINTS ARE ADDED
      NB=41
      NF=2
C...
C...  DEFINE AN INITIAL UNIFORM GRID OF NB POINTS
      CALL GRIDE(NB,XL,XU,DXB,XB)
C...
C...  THE NUMBER OF ADAPTIVE GRID POINTS, N, INITIALLY EQUALS THE
C...  NUMBER OF BASIC GRID POINTS, NB
      N=NB
C...
C...  DEFINE THE ADAPTIVE GRID INTERVAL, DXF, FROM THE BASIC GRID
C...  INTERVAL, DXB
      DXF=DXB/DFLOAT(NF)
C...
C...  DEFINE THE THRESHOLD FOR THE SECOND DERIVATIVE, ABOVE WHICH
C...  ADAPTIVE GRID POINTS ARE ADDED
      DL=20.0D+00
C...
C...  DEFINE THE MAXIMUM DISTANCES TO THE LEFT AND RIGHT OF A BASIC
C...  GRID POINT, XFL AND XFR, BEYOND WHICH ADAPTIVE GRID POINTS WILL
C...  NOT BE ADDED
      SFL=1.0D+00
      SFR=1.5D+00
      XFL=SFL*V*TP
      XFR=SFR*V*TP
C...
C...  DEFINE THE INITIAL BASIC AND ADAPTIVE GRIDS, EACH WITH NB POINTS
      DO 1 I=1,NB
C...
C...  INITIALIZE ARRAY LXB WHICH SPECIFIES THE EXISTENCE OF ADAPTIVE
C...  GRID POINTS
      LXB(I)=0
C...
C...  STORE THE SUBSCRIPTS OF THE BASIC GRID IN ARRAY NOB
      NOB(I)=I
C...
C...  DEFINE THE INITIAL CONDITIONS OF THE BASIC AND ADAPTIVE GRIDS
C...  FROM THE EXACT SOLUTION TO BURGERS EQUATION
      CB(I)=PHI(0.0D+00,XB(I))
      CA(I)=CB(I)
      XA(I)=XB(I)
1     CONTINUE
C...
C...  INITIALIZE THE DERIVATIVE CALCULATIONS IN SUBROUTINE DERV
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
C...
C...  SUBROUTINE DERV
C...
C...     (1)  COMPUTES THE FIRST AND SECOND DERIVATIVES, UX AND UXX,
C...          VIA A CALL TO SUBROUTINE NCSPLE
C...
C...     (2)  USES DERIVATIVES UX AND UXX TO ASSEMBLE BURGERS EQUATION
C...          IN THE USUAL METHOD OF LINES FORMAT
C...
C...  TYPE REAL VARIABLES AS DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
      COMMON/SYSTM1/ T0,TF,TP,H,ERROR,N
      COMMON  /T/         T
     1        /Y/   CA(150)
     2        /F/   CT(150)
     3        /B/   CB(120), CBXX(120)
     4       /SV/   XA(150)
     5       /SC/   CX(150),  CXX(150), CXXX(150)
     6        /C/        XL,        XU,         V,       VIS,
     7                   SN,        NB,        IP,        JP
      COMMON/AB1/        DL,       XFL,       XFR,        NF,
     1                  DXF,   XB(120),  LXB(120),  NOB(120)
C...
C...  APPLY THE BOUNDARY CONDITIONS FOR BURGERS EQUATION AT X = 0
C...  AND X = 1
      CA(1)=PHI(T,0.0D+00)
      CA(N)=PHI(T,1.0D+00)
C...
C...  COMPUTE THE FIRST AND SECOND SPATIAL DERIVATIVES IN BURGERS
C...  EQUATION, UX AND UXX, BY THE NORMAL CUBIC SPLINE
      CALL NCSPLE(N,XA,CA,CX,CXX,CXXX)
C...
C...  ASSEMBLE BURGERS EQUATION OVER GRID POINTS 2 TO N-1.  COMMON
C...  /SC/ ACTUALLY CONTAINS THE COEFFICIENTS OF THE CUBIC SPLINE
C...  C(X) = A0 + A1*(X - XI) + A2*(X - XI)**2 + A3*(X - XI)***3,
C...  I.E., COMMON/SC/A1(150),A2(150),A3(150), AFTER NCSPLE IS
C...  CALLED.  THUS CX(XI) = A1, CXX(XI) = 2*A2, CXXX(XI) = 6*A3,
C...  AND THE SECOND ARRAY IN /SC/ MUST BE MULTIPLIED BY 2 TO OBTAIN
C...  THE SECOND DERIVATIVE, CXX, IN BURGERS EQUATION.  ALSO, FROM
C...  THE EQUATION FOR THE SPLINE, C(XI) = A0 = CA(XI) WHERE ARRAY
C...  CA IS IN COMMON/Y/.  SUBROUTINE NCSPLE THEREFORE NUMERICALLY
C...  DIFFERENTIATES CA TO CX, CXX AND CXXX AS WELL AS EVALUATING
C...  THE COEFFICIENTS OF THE CUBIC SPLINE (THESE TWO OPERATIONS
C...  ARE EQUIVALENT)
      NM1=N-1
      DO 1 I=2,NM1
      CT(I)=2.0D+00*VIS*CXX(I)-CA(I)*CX(I)
1     CONTINUE
C...
C...  ZERO THE TIME DERIVATIVES AT I = 1, N TO INSURE THE BOUNDARY
C...  CONDITIONS ARE OBSERVED
      CT(1)=0.0D+00
      CT(N)=0.0D+00
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
C...
C...  SUBROUTINE PRINT
C...
C...     (1)  STORES THE NUMERICAL AND EXACT SOLUTIONS TO BURGERS
C...          EQUATION DURING SUCCESSIVE CALLS FOR PLOTTING
C...
C...     (2)  REDEFINES THE ADAPTIVE GRID IN ACCORDANCE WITH CHANGES
C...          IN THE NUMERICAL SOLUTION VIA A CALL TO SUBROUTINE ANUGB1
C...
C...     (3)  PLOTS THE NUMBER OF POINTS IN THE ADAPTIVE GRID VS T AT
C...          THE END OF THE RUN AND PRINTS THE AVERAGE NUMBER OF GRID
C...          POINTS DURING THE RUN
C...
C...  TYPE REAL VARIABLES AS DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
      COMMON/SYSTM1/ T0,TF,TP,H,ERROR,N
      COMMON  /T/         T
     1        /Y/   CA(150)
     2        /F/   CT(150)
     3        /B/   CB(120), CBXX(120)
     4       /SV/   XA(150)
     5       /SC/   CX(150),  CXX(150), CXXX(150)
     6        /C/        XL,        XU,         V,       VIS,
     7                   SN,        NB,        IP,        JP
      COMMON/AB1/        DL,       XFL,       XFR,        NF,
     1                  DXF,   XB(120),  LXB(120),  NOB(120)
C...
C...  DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION
      DIMENSION     GR(150),  CAA(150),    PN(51),    TN(51),
     1          CBXXP(2,51)
      REAL         IGR(150), ICAA(150),   IPN(51),   ITN(51),
     1         ICBXXP(2,51),  ICA(150),  IXA(150)
C...
C...  INCREMENT THE INDEX FOR THE PLOTTED SOLUTION
      IP=IP+1
C...
C...  COMPUTE THE RUNNING SUM OF THE NUMBER OF ODES
      SN=SN+DFLOAT(N)
C...
C...  STORE THE NUMBER OF ODES FOR SUBSEQUENT PLOTTING
      PN(IP)=DFLOAT(N)
      TN(IP)=T
      IPN(IP)=SNGL(PN(IP))
      ITN(IP)=SNGL(TN(IP))
C...
C...  CALL SUBROUTINE NCSPLE TO UPDATE THE SPLINE COEFFICIENTS USED
C...  IN THE NEW ADAPTIVE GRID
      CALL NCSPLE(N,XA,CA,CX,CXX,CXXX)
C...
C...  STORE THE DEPENDENT VARIABLE AND ITS SECOND DERIVATIVE ON THE
C...  ADAPTIVE GRID, CA(I) AND CXX(I), AT THE BASIC GRID POINTS WITH
C...  SUBSCRIPT K IN PREPARATION FOR UPDATING THE ADAPTIVE GRID.  THE
C...  FACTOR OF 2 MULTIPLYING CXX(K) IS EXPLAINED IN SUBROUTINE DERV
      DO 1 I=1,NB
      K=NOB(I)
      CB(I)=CA(K)
      CBXX(I)=2.0D+00*CXX(K)
1     CONTINUE
C...
C...  STORE THE MINIMUM AND MAXIMUM VALUES OF THE SECOND DERIVATIVE,
C...  CBXX, FOR SUBSEQUENT PLOTTING
      CBMAX=DABS(CBXX(1))
      CBMIN=DABS(CBXX(1))
      DO 10 I=2,NB
      CBMAX=DMAX1(DABS(CBXX(I)),CBMAX)
      CBMIN=DMIN1(DABS(CBXX(I)),CBMIN)
10    CONTINUE
      CBXXP(1,IP)=CBMAX
      CBXXP(2,IP)=CBMIN
      ICBXXP(1,IP)=SNGL(CBXXP(1,IP))
      ICBXXP(2,IP)=SNGL(CBXXP(2,IP))
C...
C...  STORE THE EXACT SOLUTION FOR SUBSEQUENT PLOTTING
      IF((IP-6)*(IP-16)*(IP-26)*(IP-36)*(IP-46).NE.0)GO TO 3
      JP=JP+1
      DO 2 I=1,N
      CAA(I)=PHI(T,XA(I))
      GR(I)=-0.05D+00*DFLOAT(5-JP)
      ICAA(I)=SNGL(CAA(I))
      IGR(I)=SNGL(GR(I))
      ICA(I)=SNGL(CA(I))
      IXA(I)=SNGL(XA(I))
2     CONTINUE
C...
C...  PLOT THE NUMERICAL SOLUTION WITH THE ADAPTIVE GRID SUPERIMPOSED
C...  AT SUCCESSIVE VALUES OF THE INDEPENDENT VARIABLE T.  INTEGER
C...  JP SPECIFIES THAT ONLY ADDITIONAL CURVES ARE PLOTTED AFTER THE
C...  FIRST CALL TO QIKPLT
      IM=-N
      IF(JP.GT.1)GO TO 4
C...
C...  ******************************************************************
C...  THE FOLLOWING CODING CALLS A CONTINUOUS (CALCOMP) PLOTTER.  IT
C...  IS NOT STANDARD FORTRAN, AND THEREFORE WILL PROBABLY HAVE TO BE
C...  CHANGED FOR OTHER COMPUTERS AND PLOTTERS
C...  CALL QIKSET(4.,0.,0.25,7.,-0.3,0.2)
C...  CALL QIKSAX(3,3)
C...  CALL QIKFORM(2,2)
C...  CALL QIKPLT(IXA,ICAA,N,3H*X*,8H*U(X,T)*,18H*BURGERS EQUATION*,1)
C...  CALL PLOT(-5.,1.,-3)
C...  CALL QLINE(IXA,ICA,IM,4)
C...  CALL QLINE(IXA,IGR,IM,-13)
      GO TO 6
4     CONTINUE
C...  CALL QLINE(IXA,ICAA,N,M0)
C...  CALL QLINE(IXA,ICA,IM,4)
C...  CALL QLINE(IXA,IGR,IM,-13)
C...  ******************************************************************
C...
C...  PRINT THE NUMERICAL SOLUTION CORRESPONDING TO THE PLOTTED SOLUTION
6     WRITE(NO,7)T,N,(XA(I),I=1,N)
7     FORMAT(1H ,//,5H T = ,D11.3,3X,5H N = ,I3,/,3H  X,/,(10D11.3))
      WRITE(NO,8)(CAA(I),I=1,N)
8     FORMAT(1H ,7H U ANAL,/,(10D11.3))
      WRITE(NO,9)( CA(I),I=1,N)
9     FORMAT(1H ,7H U NUM ,/,(10D11.3))
C...
C...  CONTINUE WITH SUBROUTINE PRINT, PARTICULARLY THE REDEFINITION OF
C...  THE ADAPTIVE GRID
3     CONTINUE
C...
C...  CALL SUBROUTINE ANUGB1 TO REDEFINE THE ADAPTIVE GRID
      CALL ANUGB1(NB,CB,CBXX,N,XA,CA)
C...
C...  AT THE END OF THE RUN, PLOT THE NUMBER OF ODES VS T
      IF(IP.LT.51)RETURN
      CALL TPLOTS(1,IP,TN,PN)
C...
C...  COMPUTE AND PRINT THE AVERAGE NUMBER OF ODES
      AN=SN/DFLOAT(IP)
      WRITE(NO,5)AN
5     FORMAT(1H ,//,
     1 7H N VS T,3X,8HAVG N = ,F6.2)
C...
C...  PLOT THE NUMBER OF ODES, THE MINIMUM AND MAXIMUM VALUES OF THE
C...  SECOND DERIVATIVE VS T
C...  ******************************************************************
C...  THE FOLLOWING CODING CALLS A CONTINUOUS (CALCOMP) PLOTTER.  IT
C...  IS NOT STANDARD FORTRAN, AND THEREFORE WILL PROBABLY HAVE TO BE
C...  CHANGED FOR OTHER COMPUTERS AND PLOTTERS
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(4.,0.,0.25,7.,40.,2.5)
C...  CALL QIKPLT(ITN,IPN,IP,3H*T*,3H*N*,16H*NUMBER OF ODES*,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKSET(4.,0.,0.25,7.,0.,400.)
C...  CALL QIKPLT(ITN,ICBXXP,IP,3H*T*,6H*CBXX*,19H*SECOND DERIVATIVE*,2)
C...  ******************************************************************
C...
      RETURN
      END
      DOUBLE PRECISION FUNCTION PHI(T,X)
C...
C...  FUNCTION PHI(T,X) COMPUTES THE EXACT SOLUTION TO BURGERS
C...  EQUATION.  IT IS USED IN THREE PLACES
C...
C...     (1)  TO PROVIDE THE EXACT SOLUTION IN SUBROUTINE PRINT FOR
C...          COMPARISON WITH THE NUMERICAL SOLUTION.
C...
C...     (2)  TO PROVIDE THE INITIAL CONDITION IN SUBROUTINE INITAL.
C...
C...     (3)  TO PROVIDE THE BOUNDARY CONDITIONS IN SUBROUTINE DERV.
C...
C...  ARGUMENT LIST
C...
C...      T   INITIAL-VALUE INDEPENDENT VARIABLE IN BURGERS EQUATION
C...          (INPUT)
C...
C...      X   BOUNDARY-VALUE INDEPENDENT VARIABLE IN BURGERS EQUATION
C...          (INPUT)
C...
C...  THE VALUE OF THE FUNCTION (PHI) IS THE EXACT SOLUTION TO BURGERS
C...  EQUATION AT X AND T.
C...
C...  TYPE REAL VARIABLES AS DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
      COMMON  /C/        XL,        XU,         V,       VIS
      A=(0.05D+00/VIS)*(X-0.5D+00+4.95D+00*T)
      B=(0.25D+00/VIS)*(X-0.5D+00+0.75D+00*T)
      C=( 0.5D+00/VIS)*(X-0.375D+00)
      EA=DEXP(-A)
      EB=DEXP(-B)
      EC=DEXP(-C)
      PHI=(0.1D+00*EA+0.5D+00*EB+EC)/(EA+EB+EC)
      RETURN
      END
NMOL SOLUTION OF BURGERS EQUATION BY SPLINE APPROXIMATION AND LSODE
0.        1.0       0.02
   41 1000    8    1  ABS 0.001
END OF RUNS
