      SUBROUTINE INITAL
C...
C...  KROGH'S ODE PROBLEM
C...
C...  KROGH'S PROBLEM IS A NONLINEAR ODE PROBLEM WITH AN EXACT
C...  SOLUTION.  THE NUMBER OF ODES AND THE STIFFNESS OF THESE
C...  EQUATIONS CAN BE VARIED BY THE CHOICE OF PARAMETERS.  IN
C...  THIS CASE, N = 20, AND THE STIFFNESS IS SET BY CONSTANTS
C...  READ FROM A DATA FILE.
C...
C...  THE NONLINEAR ODE PROBLEM PROPOSED BY KROGH (GEAR (1), PP.
C...  218-220), IS AN IDEAL TEST PROBLEM FOR ODE SOFTWARE FOR THE
C...  FOLLOWING REASONS:
C...
C...     (1)  THE ORDER OF THE PROBLEM (NUMBER OF ODES) CAN BE
C...          SELECTED BY THE USER (IN THE PRESENT CASE, N = 20).
C...
C...     (2)  THE STIFFNESS OF THE ODES CAN BE SELECTED BY THE USER.
C...          THIS POINT REQUIRES SOME CLARIFICATION SINCE THE ODES
C...          ARE NONLINEAR, BUT STIFFNESS IS A LINEAR CONCEPT BASED
C...          ON EIGENVALUES (A STIFF SYSTEM IS GENERALLY DEFINED AS:
C...
C...                 (PROBLEM TIME)*(LARGEST EIGENVALUE)
C...
C...                       IS MUCH GREATER THAN 1)
C...
C...          GENERALLY, THE CONCEPT OF STIFFNESS APPLIED TO NONLINEAR
C...          PROBLEMS MEANS APPLIED TO THE LINEARIZED VERSION OF THE
C...          NONLINEAR EQUATIONS, FOR WHICH EIGENVALUES CAN BE COM-
C...          PUTED.  HOWEVER, SINCE LINEARIZATION IS PERFORMED AROUND
C...          A BASE POINT, THE EIGENVALUES OF THE LINEARIZED SYSTEM
C...          WILL CHANGE WITH THE BASE POINT, I.E., THE LOCAL EIGEN-
C...          VALUES OF A NONLINEAR SYSTEM ARE NOT CONSTANT.
C...
C...          THE NONLINEAR ODES IN THE KROGH PROBLEM HAVE IN THEIR
C...          ANALYTICAL SOLUTION EXPONENTIALS OF THE FORM EXP(B(I)*T)
C...          WHERE T IS THE INDEPENDENT VARIABLE AND THE B(I), I =
C...          1, 2,..., N ARE PARAMETERS SELECTED BY THE USER.  THUS
C...          THE B(I) ARE ANALOGOUS TO EIGENVALUES, BUT THE EXPONENT-
C...          IALS (EXP(B(I)*T)) APPEAR NONLINEARLY IN THE ANALYTICAL
C...          SOLUTION (RATHER THAN AS A LINEAR COMBINATION IN THE
C...          ANALYTICAL SOLUTION OF A SYSTEM OF CONSTANT COEFFICIENT
C...          LINEAR ODES).
C...
C...          IN THIS SENSE, THE STIFFNESS OF THE KROGH TEST PROBLEM
C...          CAN BE SPECIFIED BY THE USER BY SELECTING THE B(I).
C...
C...     (3)  AS THE PRECEDING DISCUSSION INDICATES, THE KROGH TEST
C...          HAS AN ANALYTICAL (EXACT) SOLUTION, WHICH CAN BE USED
C...          AS A STANDARD IN EVALUATING A NUMERICAL INTEGRATOR, I.E.,
C...          THE ERRORS IN NUMERICAL SOLUTIONS CAN BE COMPUTED.
C...
C...  ALL OF THE FEATURES OF THE KROGH TEST PROBLEM ARE ILLUSTRATED IN
C...  THE FOLLOWING CODING FOR TWO NONSTIFF CASES:
C...
C...      (1)  WE TAKE B(I) = -I, I = 1, 2,..., 20
C...
C...          THESE CONSTANTS ARE READ FROM A DATA FILE IN SUBROUTINE
C...          INITAL AT THE BEGINNING OF THE NUMERICAL SOLUTION, AND
C...          THE INTEGRATION IS PERFORMED BY BY THE SECOND-ORDER
C...          ADAMS-BASHFORTH METHOD
C...
C...      (1)  WE TAKE B(I) = I, I = 1, 2,..., 20 (NOTE THE ONLY CHANGE
C...           IS IN THE SIGN OF B(I)).
C...
C...  THE MATHEMATICAL DETAILS OF THE PROBLEM ARE GIVEN IN REFERENCE
C...  (1) CITED PREVIOUSLY.  THEY CAN ALSO BE INFERRED FROM THE
C...  FOLLOWING CODING.
C...
      COMMON /T/        T
      COMMON /Y/    Y(20)
      COMMON /F/ DYDT(20)
      COMMON /C/ BETA(20), U(20,20), B(20,20), N
      COMMON /IO/  NI, NO
C...
C...  READ THE NUMBER OF ODES AND THE CONSTANT ARRAY, BETA, WHICH
C...  DEFINES THE PROBLEM TIME SCALE AND STIFFNESS
      READ(NI,100)N
      DO 10 J=1,N
         READ(NI,101)BETA(J)
  10  CONTINUE
      WRITE(NO,102)N,(BETA(J),J=1,N)
C...
C...  INITIAL CONDITIONS
      DO 20 J=1,N
         Y(J) = -1.0
  20  CONTINUE
C...
C...  SET THE CONSTANT ARRAY, U
      TERM = 2.0/FLOAT(N)
      DO 30 I=1,N
      DO 40 J=1,N
         U(I,J) = TERM
  40  CONTINUE
         U(I,I) = TERM - 1.0
  30  CONTINUE
C...
C...  SET THE CONSTANT ARRAY, B
      DO 50 J=1,N
      DO 60 I=1,N
         SUM = 0.0
      DO 70 K=1,N
         SUM = SUM + U(I,K)*BETA(K)*U(K,J)
  70  CONTINUE
         B(I,J) = SUM
  60  CONTINUE
  50  CONTINUE
C...
C...  FORMATS ARE LISTED HERE
 100  FORMAT(I5)
 101  FORMAT(F20.5)
 102  FORMAT(/,11X,'N = ',I3,10X,'BETA',20(/,17X,F15.3))
      RETURN
      END
      SUBROUTINE DERV
      COMMON /T/        T
      COMMON /Y/    Y(20)
      COMMON /F/ DYDT(20)
      COMMON /C/ BETA(20), U(20,20), B(20,20), N
      DIMENSION W(20)
C...
C...  CALCULATE THE ARRAY, W
      DO 10 I=1,N
         SUM = 0.0
      DO 20 J=1,N
         SUM = SUM + U(I,J)*Y(J)
  20  CONTINUE
         W(I) = SUM
  10  CONTINUE
C...
C...  CALCULATE THE N ODE DERIVATIVES
      DO 30 I=1,N
         SUM = 0.0
      DO 40 J=1,N
         SUM = SUM + U(I,J)*W(J)*W(J) - B(I,J)*Y(J)
  40  CONTINUE
         DYDT(I) = SUM
  30  CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON /T/        T
      COMMON /Y/    Y(20)
      COMMON /F/ DYDT(20)
      COMMON /C/ BETA(20), U(20,20), B(20,20), N
      DIMENSION W(20), YEX(20)
C...
C...  ABSOLUTE DIMENSIONING OF THE ARRAYS REQUIRED BY SUBROUTINE JMAP
C...  (A, X, XOLD, F, FOLD) AND SUBROUTINE EIGEN (WI, WR, Z, RW, IW)
      REAL
     1           A(20,20), WR(20), WI(20), Z(20,20), RW(20),
     2                              X(20), XOLD(20),
     3                              F(20), FOLD(20)
      INTEGER                      IW(20)
C...
C...  THE FOLLOWING EQUIVALENCE FACILITATES TRANSFER OF THE DEPENDENT
C...  VARIABLE VECTOR IN /Y/ TO ARRAY Y, AND THE DERIVATIVE VECTOR IN
C...  IN /F/ TO ARRAY F FOR USE IN THE SUBSEQUENT CALL TO SUBROUTINE
C...  JMAP
      EQUIVALENCE (Y(1),X(1)),(DYDT(1),F(1))
C...
C...  INITIALIZE A COUNTER FOR THE CALLS TO SUBROUTINES JMAP AND EIGEN
      DATA IP/0/
      IP=IP+1
C...
C...  MONITOR THE PROGRESS OF THE SOLUTION ON THE SCREEN
      WRITE(*,*)' IP = ',IP
C...
C...  MAP THE JACOBIAN MATRIX OF THE ODE SYSTEM DEFINED IN SUBROUTINE
C...  DERV, AND COMPUTE ITS TEMPORAL EIGENVALUES.  SUBROUTINE EIGEN
C...  CALLS A SERIES OF EISPACK ROUTINES TO COMPUTE THE TEMPORAL
C...  EIGENVALUES, AND OPTIONALLY, THE EIGENVECTORS OF THE 20-ODE
C...  SYSTEM JACOBIAN MATRIX (THE CALLS TO JMAP AND EIGEN ARE MADE
C...  ONLY AT THE BEGINNING AND END OF EACH RUN TO REDUCE THE VOLUME
C...  OF OUTPUT)
      IF(IP.EQ.1)THEN
         CALL  JMAP(N,A, X,XOLD,F,FOLD   )
         CALL EIGEN(N,A,WR,  WI,Z,  RW,IW)
C...
C...  PRINT A HEADING FOR THE SOLUTION
         WRITE(NO,2)
2        FORMAT(//,3X,'TIME',5X,'Y(1)',3X,'YEX(1)',5X,'PER1',
     1                       4X,'Y(20)',2X,'YEX(20)',4X,'PER20',
     2                       3X,'PERMAX')
      END IF
C...
C...  CALCULATE THE N EXACT SOLUTIONS
      DO 10 J=1,N
         BETAT = BETA(J)*T
      IF (BETAT.GT.200.) THEN
         W(J) = 0.0
      ELSE IF (BETAT.LT.-200.) THEN
         W(J) = BETA(J)
      ELSE
         W(J) = BETA(J)/(1.0 - (1.0 + BETA(J))*EXP(BETAT))
      END IF
  10  CONTINUE
      DO 20 I=1,N
         SUM = 0.0
      DO 30 J=1,N
         SUM = SUM + U(I,J)*W(J)
  30  CONTINUE
         YEX(I) = SUM
  20  CONTINUE
C...
C...  CALCULATE THE MAXIMUM ERROR
      ERRMAX = ABS((Y(1) - YEX(1))/YEX(1))*100.
      DO 40 J=2,N
         ERR = ABS((Y(J) - YEX(J))/YEX(J))*100.
         IF (ERR.GT.ERRMAX) ERRMAX = ERR
  40  CONTINUE
C...
C...  PRINT THE NUMERICAL AND EXACT SOLUTIONS FOR I = 1 AND N, AND
C...  THE CORRESPONDING ERRORS, AND THE MAXIMUM ERROR
      PER1=ABS((Y(1)-YEX(1))/YEX(1))*100.
      PERN=ABS((Y(N)-YEX(N))/YEX(N))*100.
      PERMAX=ERRMAX
      WRITE(NO,3)T,Y(1),YEX(1),PER1,Y(N),YEX(N),PERN,PERMAX
3     FORMAT(F7.2,7F9.4)
C...
C...  MAP THE JACOBIAN MATRIX OF THE ODE SYSTEM DEFINED IN SUBROUTINE
C...  DERV, AND COMPUTE ITS TEMPORAL EIGENVALUES AT THE END OF THE RUN
      IF(IP.EQ.11)THEN
         CALL  JMAP(N,A, X,XOLD,F,FOLD   )
         CALL EIGEN(N,A,WR,  WI,Z,  RW,IW)
      END IF
C...
C...  RESET IP AT THE END OF EACH RUN
      IF(IP.EQ.11)IP=0
      RETURN
      END





