      SUBROUTINE INITAL
C...
C...     DYI/DT + I*YI = YI-1, YI(0) = 0, (Y0 = 1)           (1)(2)
C...
      PARAMETER (N=5)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(N)      
     2      /F/   DYDT(N)
C...
C...  INITIAL CONDITION (2)
      DO 1 I=1,N
      Y(I)=0.
1     CONTINUE
C...
C...  COMPUTE THE INITIAL DERIVATIVES
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
      PARAMETER (N=5)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(N)      
     2      /F/   DYDT(N)
C...
C...  EQUATION (1)
      DO 1 I=1,N
      IF(I.EQ.1)THEN
         DYDT(I)=-FLOAT(I)*Y(I)+1.0
      ELSE
         DYDT(I)=-FLOAT(I)*Y(I)+Y(I-1)
      END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      PARAMETER (N=5)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(N)      
     2      /F/   DYDT(N)
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(5,5), WR(5), WI(5), Z(5,5),   RW(5),
     2                                   X(5), XOLD(5),
     3                                   F(5), FOLD(5)
      INTEGER                   IW(5)
C...
C...  ARRAY FOR THE COEFFICIENTS IN THE EXACT SOLUTION
      REAL      AC(0:5)
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 AN INTEGER COUNTER TO CONTROL THE PRINTING
      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 5-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)
         WRITE(NO,5)
5        FORMAT(1H ,//)
C...
C...  PRINT A HEADING FOR THE OUTPUT
         WRITE(NO,2)
2        FORMAT(1H ,//,9X,'T'10X,'Y5(T)',5X,'ANAL Y5(T)',11X,'DIFF')
C...
C...  COMPUTE THE ANALYTICAL SOLUTION AND THE PERCENT DIFFERENCE WITH
C...  THE NUMERICAL SOLUTION.  THE COEFFIENTS IN THE ANALYTICAL SOLU-
C...  TION ARE CALCULATED ONLY ONCE
         DO 3 I=0,N
            AC(I)=1.0
         DO 3 J=0,N
            IF(I.NE.J)THEN
            AC(I)=AC(I)*1.0/FLOAT(J-I)
            END IF
3     CONTINUE
      END IF
C...
C...  SUM THE EXPONENTIALS IN THE EXACT SOLUTION
      YA=0.
      DO 4 I=0,N
      YA=YA+AC(I)*EXP(-FLOAT(I)*T)
4     CONTINUE
C...
C...  PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS
      IF(IP.EQ.1)THEN
         DIFF=0.
      ELSE
         DIFF=(Y(N)-YA)/YA*100.
      END IF
      WRITE(NO,1)T,Y(N),YA,DIFF
1     FORMAT(F10.3,3E15.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
C...  RUN
      IF(IP.EQ.11)THEN
         CALL  JMAP(N,A, X,XOLD,F,FOLD   )
         CALL EIGEN(N,A,WR,  WI,Z,  RW,IW)
         WRITE(NO,5)
      END IF
C...
C...  RESET THE COUNTER IP AT THE END OF EACH RUN
      IF(IP.EQ.11)IP=0
      RETURN
      END

