*DECK REACT35
      SUBROUTINE INITAL
C...
C...  SOLUTION OF A SECOND ORDER, NONLINEAR, TWO POINT BOUNDARY VALUE
C...  ORDINARY DIFFERENTIAL EQUATION BY QUASILINEARIZATION AND THE
C...  METHOD OF FALSE TRANSIENTS - APPLICATION IN CHEMICAL KINETICS
C...
C...  CONSIDER THE HEATING OF A SOLID IN WHICH AN EXOTHERMIC REACTION
C...  TAKES PLACE.  A STEADY STATE ENERGY BALANCE ON A SECTION OF THE
C...  SOLID OF THICKNESS OF DX GIVES:
C...
C...     A*Q(X) - A*Q(X+DX) = HR*A*DX*RHO*K0*EXP(-E/(R*T))
C...
C...  WHERE
C...
C...     Q       HEAT FLUX BY CONDUCTION (CAL/CM**2-SEC)
C...
C...     T       ABSOLUTE TEMPERATURE (K)
C...
C...     X       DISTANCE ALONG THE SOLID (CM)
C...
C...     A       CROSS SECTIONAL AREA OF THE SECTION (CM**2)
C...
C...     HR      HEAT OF REACTION (CAL/GM)
C...
C...     RHO     DENSITY (GM/CM**3)
C...
C...     K0      PRE-EXPONENTIAL FACTOR (1/SEC)
C...
C...     E       ACTIVATION ENERGY (CAL/GM MOL)
C...
C...     R       GAS CONSTANT (CAL/GM MOL-K)
C...
C...  DIVISION BY A*DX, FOLLOWED BY DX ---> 0 GIVES
C...
C...     -Q  = HR*RHO*K0*EXP(-E/*R*T))                          (1)
C...       X
C...
C...  WHERE A LETTER SUBSCRIPT DENOTES A DERIVATIVE, E.G., Q  = DQ/DX.
C...                                                        X
C...
C...  EQUATION (1) CONTAINS TWO DEPENDENT VARIABLES, Q AND T, SO THAT
C...  A SECOND EQUATION IS REQUIRED
C...
C...     Q = -K*T                                               (2)
C...             X
C...
C...  EQUATION (2) IS FOURIER'S FIRST LAW.  SUBSTUTION OF EQUATION (2)
C...  IN (1) GIVES
C...
C...     T   = (HR*RHO/K)*K0*EXP(-E/(R*T))                      (3)
C...      XX
C...
C...  EQUATION (3) IS A NONLINEAR, BOUNDARY VALUE ORDINARY DIFFERENTIAL
C...  EQUATION, WHICH WILL BE INTEGRATED NUMERICALLY SUBJECT TO THE TWO
C...  BOUNDARY CONDITIONS
C...
C...     T(0) = TA, T (L) = 0                                (4)(5)
C...                 X
C...
C...     TA      AMBIENT TEMPERATURE (K)
C...
C...     L       LENGTH OF THE SOLID (CM)
C...
C...  IF THE SECOND DERIVATIVE IN EQUATION (3) IS REPLACED WITH A SECOND
C...  ORDER APPROXIMATION
C...
C...     T   = (T(I+1) - 2*T(I) + T(I-1))/DX**2
C...      XX
C...
C...  EQUATION (3) BECOMES
C...
C...     (T(I+1) - 2*T(I) + T(I-1))/DX**2 - D*EXP(-E/(R*T(I)) = 0    (6)
C...
C...  WHERE
C...
C...     D       HR*RHO*K0/K (K/CM**2)
C...
C...     DX      GRID SPACING (CM)
C...
C...  IF EQUATION (6) IS NOW WRITTEN FOR THE N+1 GRID POINTS, I = 0,
C...  1, 2,..., N, A SYSTEM OF NONLINEAR TRIDIAGONAL TRANSCENDENTAL
C...  EQUATIONS RESULTS, WHICH CAN BE APPROXIMATED BY LINEAR FUNCTIONS
C...  AS EXPLAINED IN CLASS.
C...
C...  COMPUTE A SOLUTION TO THE NONLINEAR PROBLEM, I.E., T(X), USING
C...  SUBROUTINES DECOMP AND SOLVE, AND TRIDAG, FOR:
C...
C...     HR = -10000 CAL/GM, RHO = 1 GM/CM**3, K0 = 30.0 (1/SEC)
C...
C...     K = 0.01 CAL-CM/SEC-CM**2-K, E = 10000 CAL/GM MOL, R = 1.987
C...
C...     CAL/GM MOL-K, TA = 350 K, L = 1 CM
C...
C...  (6) IS A SET OF NONLINEAR, TRIDIAGONAL TRANSCENDENTAL EQUATIONS.
C...  WE WILL COMPUTE A SOLUTION BY NEWTON'S METHOD, WHICH IS ALSO
C...  CALLED QUASILINEARIZATION.  IF THE NONLINEAR EXPONENTIAL TERM
C...  IS EXPANDED IN A TAYLOR SERIES
C...
C...                   K+1                  K
C...     EXP(-E/(R*T(I)   ) = EXP(-E/(R*T(I) )) +
C...
C...                                          K            K+1       K
C...                        + D(EXP(-E/(R*T(I) )))/DT*(T(I)    - T(I) )
C...
C...                                          K+1
C...                        = A(I) + B(I)*T(I)    + C(I)
C...
C...  WHERE
C...
C...     K       ITERATION COUNTER
C...
C...                           K
C...     A(I)    EXP(-E/(R*T(I) ))
C...
C...                             K
C...     B(I)    D(EXP(-E/(R*T(I) )))/DT =
C...
C...                           K
C...             EXP(-E/(R*T(I) ))*(E/(R*T(I)**2))
C...
C...                        K
C...     C(I)    B(I)*(-T(I) )
C...
C...  THEN EQUATIONS (6) CAN BE WRITTEN IN APPROXIMATE LINEAR FORM AS
C...  (WITH THE ITERATION COUNTER K DROPPED)
C...
C...     - (2 + (DX**2)*D*B(1))*T(1) + T(2) =
C...
C...          - TA + (DX**2)*D*(A(1) + C(1))
C...
C...     T(1) - (2 + (DX**2)*D*B(2))*T(2) + T(3) =
C...
C...                (DX**2)*D*(A(2) + C(2))
C...       .                                      .
C...       .                                      .             (7)
C...       .                                      .
C...
C...     T(N-2) - (2 + (DX**2)*D*B(N-1))*T(N-1) + T(N) =
C...
C...                (DX**2)*D*(A(N-1) + C(N-1))
C...
C...     2*T(N-1) - (2 + (DX**2)*D*B(N))*T(N) =
C...
C...                (DX**2)*D*(A(N) + C(N))
C...
C...  NOTE THAT A, B AND C IN EACH EQUATION ARE FUNCTIONS OF I, AND
C...  THEREFORE MUST BE EVALUATED FOR EACH EQUATION.
C...
C...  (7) IS A SET OF LINEAR, TRIDIAGONAL ALGEBRAIC EQUATIONS WHICH
C...  CAN BE SOLVED BY A BANDED SOLVER.  IT WILL BE NECESSARY TO ASSUME
C...  A TRIAL SOLUTION FOR ITERATION K = 0.  THEN (7) CAN BE SOLVED
C...  FOR THE K = 1 SOLUTION, AND THIS SOLUTION CAN BE USED IN PLACE OF
C...  THE ASSUMED K = 0 SOLUTION TO COMPUTE THE K = 2 SOLUTION, ETC.,
C...  UNTIL CONVERGENCE IS ACHIEVED.
C...
C...  TWO SOLUTIONS TO EQUATION (3) ARE COMPUTED BY THIS METHOD FOR
C...  N = 10 AND N = 20.  SOLUTIONS ARE ALSO COMPUTED BY SUBROUTINE
C...  SNSQE WHICH CAN ACCEPT THE NONLINEAR EQUATIONS (EQUATIONS (6))
C...  DIRECTLY.
C...
C...  A SECOND APPROACH TO THE SOLUTION OF EQUATION (3) IS TO CONVERT
C...  IT TO A TRANSIENT (TIME DEPENDENT) EQUATION BY ADDING A TIME
C...  DERIVATIVE
C...
C...     RHO*CP*T   = K*T   - (HR*RHO)*K0*EXP(-E/(R*T))       (8)
C...             TM      XX
C...
C...  WHICH CAN BE INTEGRATED UNTIL THE TIME DERIVATIVE, T  ---+ 0.
C...                                                      TM
C...  THIS APPROACH, WHICH IS ALSO IMPLEMENTED IN THE FOLLOWING
C...  CODING, IS CALLED THE METHOD OF FALSE TRANSIENTS.
C...
C...  SUBROUTINES DECOMP AND SOLVE, TRIDAG, AND SNSQE WERE TAKEN FROM
C...  THE FOLLOWING REFERENCES:
C...
C...  (1)  DECOMP AND SOLVE: FORSYTHE, GEORGE E., MICHAEL A. MALCOLM,
C...       AND CLEVE B. MOLER, COMPUTER METHODS FOR MATHEMATICAL COMPU-
C...       TATIONS, PRENTICE-HALL, ENGLEWOOD CLIFFS, NJ, 1977, PP48-55.
C...
C...  (2)  TRIDAG: PRESS, WILLIAM H., BRIAN P. FLANNERY, SAUL A. 
C...       TEUKOLSKY AND WILLIAM T. VETTERLING, NUMERICAL RECIPES,
C...       CAMBRIDGE UNIVERSITY PRESS, CAMBRIDGE, 1986, P40
C...
C...  (3)  KANAHER, DAVID, CLEVE MOLER AND STEPHEN NASH, NUMERICAL
C...       METHODS AND SOFTWARE, PRENTICE HALL, ENGLEWOOD CLIFFS,
C...       1989, PP258-271
C...
C...  OTHER REFERENCES FOR THE NUMERICAL SOLUTION OF TWO POINT, BOUND-
C...  ARY VALUE ODES INCLUDE:
C...
C...  (4)  FOX, L., THE NUMERICAL SOLUTION OF TWO-POINT BOUNDARY
C...       PROBLEMS IN ORDINARY DIFFERENTIAL EQUATIONS, DOVER, NY,
C...       1990
C...
C...  (5)  KELLER, HERBERT, NUMERICAL METHODS FOR TWO-POINT BOUNDARY-
C...       VALUE PROBLEMS, BLAISDELL, WALTHAM, MA, 1968
C...
C...  (6)  LEE, E. STANLEY, QUASILINEARIZATION AND INVARIANT IMBEDDING,
C...       ACADEMIC PRESS, NEW YORK, 1968
C...
C...  TWO POINT, BOUNDARY VALUE PROBLEMS CAN ALSO BE SOLVED NUMERICALLY
C...  BY THE QUALITY, FORTRAN ROUTINE COLSYS:
C...
C...  (7)  ASCHER, U., J. CHRISTIANSEN AND R. D. RUSSELL, COLLOCATION
C...       SOFTWARE FOR BOUNDARY VALUE ODES, ACM TRANSACTIONS ON MATHE-
C...       MATICAL SOFTWARE, JUNE, 1981
C...
C...  INPUT/OUTPUT (I/O) UNIT NUMBERS ARE IN COMMON/IO/
      COMMON/IO/ NI, NO
C...
C...  DSS/2 COMMON AREA
      COMMON/T/   TIME,  NSTOP,  NORUN
     1      /Y/ TC(21)
     2      /F/ TT(21)
     3      /S/ TX(21),TXX(21)
     4      /C/      N,     IP
C...
C...  COMMON AREA FOR ACCESS IN THE FUNCTION CALLED BY SUBROUTINE 
C...  SNSQE, AND THE USUAL DSS/2 SUBROUTINES (INITAL, DERV, PRINT)
      COMMON/C1/ T0(20), D, E, R, DXS, TA, DF
C...
C...  DIMENSION THE ARRAYS REQUIRED BY SUBROUTINES DECOMP AND SOLVE
      DIMENSION AFSM(20,20), BRHS(20), WORK(20)
      INTEGER      IPVT(20)
C...
C...  DIMENSION THE ARRAYS REQUIRED BY SUBROUTINE TRIDAG
      DIMENSION AL(20), BM(20), CU(20)
C...
C...  DIMENSION THE ARRAYS REQUIRED BY SUBROUTINE SNSQE
      DIMENSION T(20), FVEC(20), WA(750)
C...
C...  EXTERNAL THE ROUTINE CALLED BY SNSQE TO EVALUATE THE FUNCTIONS.
C...  JAC TO EVALUATE THE JACOBIAN MATRIX IS NOT CALLED BY SNSQE SINCE
C...  IOPT=2 IS AN INPUT TO SNSQE (THE JACOBIAN MATRIX IS DONE BY
C...  FINITE DIFFERENCES RATHER THAN ANALYTICALLY).  A DUMMY ROUTINE
C...  JAC IS PROVIDED TO SATISFY THE LOADER/LINKER
      EXTERNAL FCN, JAC
C...
C...  REAL SELECTED VARIABLES
      REAL K0, K, L
C...
C...  SET THE MODEL PARAMETERS
      HR=-10000.
      RHO=1.0
      CP=1.0
      K0=30.0
      K=0.01
      E=10000.
      R=1.987
      TA=400.
      L=1.
C...
C...  STEP THROUGH TWO CASES
      DO 7 NCASE=1,2
C...
C...  FOR THE FIRST CASE, N = 10
      IF(NCASE.EQ.1)N=10
C...
C...  FOR THE SECOND CASE, N = 20
      IF(NCASE.EQ.2)N=20
C...
C...  ****************************************************************
C...
C...  COMPUTE A SOLUTION BY SUBROUTINES DECOMP AND SOLVE USING QUASI-
C...  LINEARIZATION
C...
C...  CONSTANT D IN EQUATIONS (7)
      D=HR*RHO*K0/K
C...
C...  GRID SPACING AND ITS SQUARE
      DX=L/FLOAT(N)
      DXS=DX*DX
C...
C...  DEFINE THE INITIAL TRIAL SOLUTION, SPATIAL GRID AND INDEX
      DO 6 I=1,N
         T0(I)=400.0E0
6     CONTINUE
C...
C...  STEP THROUGH A SERIES OF ITERATIONS
      DO 8 ITER=1,100
C...
C...  SET UP THE COEFFICIENT MATRIX IN FULL STORAGE MODE (STORED IN
C...  AFSM)                            -    -       -
C...
C...  COEFFICIENT MATRIX
      DO 1 I=1,N
      DO 1 J=1,N
         AFSM(I,J)=0.0E0
         IF(J.EQ.I  )AFSM(I,J)=-(2.0E0+DXS*D*B(I))
         IF(J.EQ.I-1)AFSM(I,J)=1.0E0
         IF(J.EQ.I+1)AFSM(I,J)=1.0E0
1     CONTINUE
      AFSM(N,N-1)=2.0E0
C...
C...  RIGHT HAND SIDE VECTOR
      DO 4 I=1,N
         IF(I.EQ.1)BRHS(I)=-TA+DXS*D*(A(I)+C(I))
         IF(I.NE.1)BRHS(I)=   +DXS*D*(A(I)+C(I))
4     CONTINUE
C...
C...  PRINT THE COEFFICIENT MATRIX, RHS VECTOR (10 X 10 CASE ONLY)
      IF((NCASE.EQ.1).AND.(ITER.EQ.1))THEN
         DO 13 I=1,N
            WRITE(NO,3)(AFSM(I,J),J=1,N),BRHS(I)
3           FORMAT(//,(11F7.1))
13       CONTINUE
      END IF
C...
C...  SOLVE THE LINEAR ALGEBRAIC EQUATIONS BY SUBROUTINES DECOMP AND
C...  SOLVE
      NDIM=20
      CALL DECOMP(NDIM,N,AFSM,COND,IPVT,WORK)
C...
C...  CHECK FOR THE CONDITION OF THE COEFFICIENT MATRIX BEFORE DOING
C...  THE BACK SUBSTITUTION
      IF(COND.GT.1.0E20)THEN
         WRITE(NO,2)COND
2        FORMAT(//,' CONDITION NUMBER FROM DECOMP IS',E12.3)
         STOP
      ELSE
         CALL SOLVE(NDIM,N,AFSM,BRHS,IPVT)
      END IF
C...
C...  CHECK FOR CONVERGENCE
      DO 9 I=1,N
         IF(ABS(BRHS(I)-T0(I)).GT.0.001)GO TO 10
9     CONTINUE
C...
C...  CONVERGENCE HAS BEEN ACHIEVED, SO PRINT THE SOLUTION VECTOR AND
C...  GO ON TO THE SOLUTION BY TRIDAG
      WRITE(NO,5)N,ITER,(BRHS(I),I=1,N)
5     FORMAT(1H ,//,' SOLUTION FROM DECOMP AND SOLVE',' N = ',I5,3X,
     +              ' ITER = ',I5,//,11X,'T(X)',/,(F15.4))
      GO TO 21
C.
C...  CONVERGENCE HAS NOT BEEN ACHIEVED, SO CONTINUE THE ITERATIONS
10    DO 11 I=1,N
      T0(I)=BRHS(I)
11    CONTINUE
C...
C...  TERMINATE EXECUTION BECAUSE OF EXCESSIVE ITERATIONS
      IF(ITER.EQ.100)THEN
         WRITE(NO,12)ITER
12       FORMAT(' NUMBER OF ITERATIONS FOR DECOMP AND SOLVE = ',I4)
         STOP
      END IF
C...
C...  GO ON TO THE NEXT ITERATION
8     CONTINUE
C...
C...  ****************************************************************
C...
C...  COMPUTE A SOLUTION BY TRIDIAGONAL ROUTINE TRIDAG
C...
C...  DEFINE THE INITIAL TRIAL SOLUTION, SPATIAL GRID AND INDEX
21    DO 26 I=1,N
         T0(I)=400.0E0
26    CONTINUE
C...
C...  STEP THROUGH A SERIES OF ITERATIONS
      DO 30 ITER=1,100
C...
C...  SET UP THE COEFFICIENT MATRIX IN BAND STORAGE MODE (STORED IN
C...  ABSM)                            -    -       -
C...
C...  LOWER DIAGONAL
      DO 23 I=1,N
         IF(I.EQ.1)AL(1)=0.0E0
         IF((I.GT.1).AND.(I.LT.N))AL(I)=1.0E0
         IF(I.EQ.N)AL(N)=2.0E0
23    CONTINUE
C...
C...  MAIN DIAGONAL
      DO 24 I=1,N
         BM(I)=-(2.0+DXS*D*B(I))
24    CONTINUE
C...
C...  UPPER DIAGONAL
      DO 25 I=1,N
         IF(I.NE.N)CU(I)=1.0E0
         IF(I.EQ.N)CU(N)=0.0E0
25    CONTINUE
C...
C...  RIGHT HAND SIDE VECTOR
      DO 28 I=1,N
         IF(I.EQ.1)BRHS(1)=-TA+DXS*D*(A(I)+C(I))
         IF(I.NE.1)BRHS(I)=   +DXS*D*(A(I)+C(I))
28    CONTINUE
C...
C...  PRINT THE COEFFICIENT MATRIX, RHS VECTOR (10 X 10 CASE ONLY)
      IF((NCASE.EQ.1).AND.(ITER.EQ.1))THEN
         DO 27 I=1,N
            WRITE(NO,3)AL(I),BM(I),CU(I),BRHS(I)
27       CONTINUE
      END IF
C...
C...  SOLVE THE LINEAR ALGEBRAIC EQUATIONS BY SUBROUTINE TRIDAG
      CALL TRIDAG(AL,BM,CU,BRHS,T,N)
C...
C...  CHECK FOR CONVERGENCE
      DO 29 I=1,N
         IF(ABS(T(I)-T0(I)).GT.0.001)GO TO 32
29    CONTINUE
C...
C...  CONVERGENCE HAS BEEN ACHIEVED, SO PRINT THE SOLUTION VECTOR AND
C...  GO ON TO THE SOLUTION BY SNSQE
      WRITE(NO,22)N,ITER,(T(I),I=1,N)
22    FORMAT(1H ,//,' SOLUTION FROM TRIDAG',' N = ',I5,3X,
     +              ' ITER = ',I5,//,11X,'T(X)',/,(F15.4))
      GO TO 31
C.
C...  CONVERGENCE HAS NOT BEEN ACHIEVED, SO CONTINUE THE ITERATIONS
32    DO 33 I=1,N
        T0(I)=T(I)
33    CONTINUE
C...
C...  TERMINATE EXECUTION BECAUSE OF EXCESSIVE ITERATIONS
      IF(ITER.EQ.100)THEN
         WRITE(NO,34)ITER
34       FORMAT(' NUMBER OF ITERATIONS FOR TRIDAG = ',I4)
         STOP
      END IF
C...
C...  GO ON TO THE NEXT ITERATION
30    CONTINUE
C...
C...  ****************************************************************
C...
C...  COMPUTE A SOLUTION BY SUBROUTINE SNSQE
C...
C...  DEFINE THE PARAMETERS FOR SNSQE
C...
C...     ERROR TOLERANCE
31       TOL=1.0E-04
C...
C...     FINITE DIFFERENCE APPROXIMATION OF THE JACOBIAN
         IOPT=2
C...
C...     NO PRINTOUT OF INTERMEDIATE ITERATES
         NPRINT=0
C...
C...     LENGTH OF THE WORK ARRAY.  LWA = (3*N**2+13*N))/2 =
C...     (3*20**2+13*20))/2 = (3*400+260))/2 = 1460/2 = 730.  WE
C...     USE OVER DIMENSIONING HERE
         LWA=750
C...
C...  SET THE INITIAL ESTIMATE OF THE SOLUTION
         DO 100 I=1,N
            T(I)=400.0E0
100      CONTINUE
C...
C...  CALL SUBROUTINE SNSQE TO COMPUTE A SOLUTION TO EQUATION (6)
      CALL SNSQE(FCN,JAC,IOPT,N,T,FVEC,TOL,NPRINT,INFO,WA,LWA)
C...
C...  COMPUTE THE SUM OF SQUARES OF THE RESIDUALS
         SUM=0.0E0
         DO 103 I=1,N
            SUM=SUM+FVEC(I)**2
103      CONTINUE
C...
C...  PRINT THE SOLUTION
         WRITE(NO,101)N,SUM,INFO,(T(I),I=1,N)
101      FORMAT(1H ,//,' SOLUTION FROM SNSQE'
     1    ' N = ',I5,3X,' SUM = ',E11.3,3X,' INFO = ',I5,//,
     2    11X,'T(X)',/,(F15.4))
C...
C...  GO ON TO THE NEXT CASE
7     CONTINUE
C...
C...  ****************************************************************
C...
C...  COMPUTE A SOLUTION BY THE METHOD OF FALSE TRANSIENTS
C...
C...  INITIAL CONDITIONS
      N=21
      TC(1)=TA
      DO 102 I=2,N
C...
C...     CASE 1 (NORUN = 1), UNIFORM TEMPERATURE.  THIS INITIAL CONDI-
C...     TION GIVES A FINAL (STEADY STATE SOLUTION) WHICH IS SIGNIFI-
C...     CANTLY DIFFERENT THAN THE SOLUTION BY QUASILINEARIZATION
         IF(NORUN.EQ.1)TC(I)=400.
C...
C...     CASE 2 (NORUN = 2), TEMPERATURES FROM QUASILINEARIZATION.  THIS
C...     SOLUTION DOES NOT REMAIN AT THE QUASILINEARIZATION SOLUTION,
C...     BUT RATHER, CONVERGES TO THE SOLUTION OF CASE 1, INDICATING
C...     AN APPARENT INSTABILITY IN THE QUASILINEARIZATION SOLUTION.
C...
C...     FINALLY, IF THE INITIAL ESTIMATE OF THE SOLUTION IS CHANGED
C...     TO 350 (RATHER THAN 400), ALL FOUR APPROACHES (DECOMP AND 
C...     SOLVE, TRIDAG, SNSQE AND THE METHOD OF LINES) CONVERGE TO THE 
C...     SAME SOLUTION (TO WITHIN 1 PART IN 1.0E+05).  THIS CAN EASILY 
C...     BE VERIFIED BY USING AN EDITOR TO CHANGE =400 TO =400 IN THE 
C...     FOUR PRECEDING
C...     STATEMENTS
C...
C...                            T0(I)=400.
C...                            (FOR DECOMP AND SOLVE)
C...
C...                            T0(I)=400.
C...                            (FOR TRIDAG)
C...
C...                            T(I)=400.
C...                            (FOR SNSQE)
C...
C...                            IF(NORUN.EQ.1)TC(I)=400.
C...                            (FOR THE METHOD OF LINES)
C...
C...     THIS RESULT (THE SAME SOLUTIONS FOR AN INITIAL ESTIMATE OF
C...     350 AND DIFFERENT SOLUTIONS FOR AN INITIAL ESTIMATE OF 400,
C...     WITH THE LATTER PRODUCING AN UNSTABLE SOLUTION BY DECOMP AND
C...     SOLVE, TRIDAG AND SNSQE) GIVES SOME INDICATION OF THE INTER-
C...     ESTING AND FREQUENTLY UNPREDICTABLE PROPERTIES OF NONLINEAR 
C...     SYSTEMS.  ON THE OTHER HAND, WELL-CONDITIONED LINEAR PROBLEMS 
C...     GENERALLY HAVE UNIQUE AND PREDICTABLE SOLUTIONS.
C...
         IF(NORUN.EQ.2)TC(I)=T(I-1)
102   CONTINUE
C...
C...  RECOMPUTE CONSTANT D TO INCLUDE CP
      D=D*K/(RHO*CP)
C...
C...  THERMAL DIFFUSIVITY
      DF=K/(RHO*CP)
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      END
      FUNCTION A(I)
C...
C...  FUNCTION A COMPUTES THE FIRST TERM IN THE LINEAR APPROXIMATION
C...  OF EXP(-E/(R*T)) AS DISCUSSED PREVIOUSLY
C...
      COMMON/C1/ T0(20), D, E, R, DXS, TA, DF
      A=EXP(-E/(R*T0(I)))
      RETURN
      END
      FUNCTION B(I)
C...
C...  FUNCTION B COMPUTES THE SECOND TERM IN THE LINEAR APPROXIMATION
C...  OF EXP(-E/(R*T)) AS DISCUSSED PREVIOUSLY
C...
      COMMON/C1/ T0(20), D, E, R, DXS, TA, DF
      B=A(I)*(E/(R*(T0(I)**2)))
      RETURN
      END
      FUNCTION C(I)
C...
C...  FUNCTION C COMPUTES THE THIRD TERM IN THE LINEAR APPROXIMATION
C...  OF EXP(-E/(R*T)) AS DISCUSSED PREVIOUSLY
C...
      COMMON/C1/ T0(20), D, E, R, DXS, TA, DF
      C=B(I)*(-T0(I))
      RETURN
      END
      SUBROUTINE FCN(N,T,F,IFLAG)
C...
C...  SUBROUTINE IMPLEMENTS THE NONLINEAR EQUATIONS
C...
      COMMON/C1/ T0(20), D, E, R, DXS, TA, DF
C...
C...  DIMENSION THE ARRAYS USED BY FCN
      DIMENSION T(N), F(N)
C...
C...  IMPLEMENT THE NONLINEAR EQUATIONS
      DO 1 I=1,N
      IF(I.EQ.1)THEN
         F(1)=(T(2)-2.0E0*T(1)+TA)/DXS-D*EXP(-E/(R*T(1)))
      ELSE IF(I.EQ.N)THEN
         F(N)=(2.0E0*T(N-1)-2.0E0*T(N))/DXS-D*EXP(-E/(R*T(N)))
      ELSE
         F(I)=(T(I+1)-2.0E0*T(I)+T(I-1))/DXS-D*EXP(-E/(R*T(I)))
      END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/   TIME,  NSTOP,  NORUN
     1      /Y/ TC(21)
     2      /F/ TT(21)
     3      /S/ TX(21),TXX(21)
     4      /C/      N,     IP
      COMMON/C1/ T0(20), D, E, R, DXS, TA, DF
      REAL K0, K, L
C...
C...  BOUNDARY CONDITIONS
      TC(1)=TA
      TX(N)=0.0E0
C...
C...  SECOND ORDER SPATIAL DERIVATIVE
      CALL DSS042(0.0E0,1.0E0,N,TC,TX,TXX,1,2)
C...
C...  PDE
      DO 1 I=1,N
      TT(I)=DF*TXX(I)-D*EXP(-E/(R*TC(I)))
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/   TIME,  NSTOP,  NORUN
     1      /Y/ TC(21)
     2      /F/ TT(21)
     3      /S/ TX(21),TXX(21)
     4      /C/      N,     IP
      COMMON/C1/ T0(20), D, E, R, DXS, TA, DF
C...
C...  REAL SELECTED VARIABLES
      REAL K0, K, L
C...
C...  DIMENSION THE ARRAY FOR THE ARRHENIUS TEMPERATURE DEPENDENCY
      DIMENSION EX(21)
C...
C...  PRINT THE SOLUTION, INCLUDING THE PARTIAL DERIVATIVES AND THE
C...  ARRHENIUS TEMPERATURE DEPENDENCY
      DO 3 I=1,N
      EX(I)=-D*EXP(-E/(R*TC(I)))
3     CONTINUE
         WRITE(NO,2)TIME,( TC(I),I=1,N,5),( TT(I),I=1,N,5),
     1                   (TXX(I),I=1,N,5),( EX(I),I=1,N,5)
2        FORMAT(1H ,//,' T = ',F10.2,/,
     1          '    TC(X,T)',5F10.2,/,
     2          '    TT(X,T)',5E10.2,/,
     3          '   TXX(X,T)',5F10.2,/,
     4          '    EX(X,T)',5F10.4,//)
      RETURN
      END
      SUBROUTINE JAC
C...
C...  DUMMY ROUTINE TO SATISFY THE LOADER/LINKER (JAC IS AN EXTERNAL
C...  FOR SUBROUTINE SNSQE, BUT IT IS NOT CALLED SINCE IOPT=2 IN CALL-
C...  ING SNSQE)
      RETURN
      END

