      SUBROUTINE INITAL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/       T,   NSTOP,   NORUN
C...
C...  LINEAR PDE
      IF((NORUN.EQ.1).OR.(NORUN.EQ.2))THEN
         CALL INIT1
C...
C...  NONLINEAR PDE
      ELSE 
     +IF((NORUN.EQ.3).OR.(NORUN.EQ.4))THEN
         CALL INIT2
      END IF
      RETURN
      END      
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/       T,   NSTOP,   NORUN
C...
C...  LINEAR PDE
      IF((NORUN.EQ.1).OR.(NORUN.EQ.2))THEN
         CALL DERV1
C...
C...  NONLINEAR PDE
      ELSE 
     +IF((NORUN.EQ.3).OR.(NORUN.EQ.4))THEN
         CALL DERV2
      END IF
      RETURN
      END      
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/       T,   NSTOP,   NORUN
C...
C...  LINEAR PDE
      IF((NORUN.EQ.1).OR.(NORUN.EQ.2))THEN
         CALL PRINT1(NI,NO)
C...
C...  NONLINEAR PDE
      ELSE 
     +IF((NORUN.EQ.3).OR.(NORUN.EQ.4))THEN
         CALL PRINT2(NI,NO)
      END IF
      RETURN
      END      
      SUBROUTINE INIT1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/       T,   NSTOP,   NORUN
     1      /Y/  U1(NX),  U2(NX)
     2      /F/ U1T(NX), U2T(NX)
     3      /S/ U1X(NX),U1XX(NX),   X(NX)
     4      /C/       A,       B,       C,       D,
     5               C1,      C2,      C3
     6      /I/      IP
C...
C...  PROBLEM PARAMETERS
      A=1.0D0
      B=1.0D0
      C=1.0D0
C...
C...  MAXWELL'S EQUATION
      IF(NORUN.EQ.1)THEN
         D=0.0D0
C...
C...  TELEGRAPH EQUATION 
      ELSE IF(NORUN.EQ.2)THEN
         D=-1.0D0
      END IF
C...
C...  PRECOMPUTE THE CONSTANTS IN THE ANALYTICAL SOLUTION FOR USE IN
C...  FUNCTION UANAL
      C3=4.0D0*DATAN(1.0D0)
      C1=-B/(2.0D0*A)
      C2=DSQRT((1.0D0/A)*((B**2)/(4.0D0*A)
     1                   -(B**2)/(2.0D0*A)
     2                   +C*(C3**2)-D))
C...
C...  SPATIAL GRID
      DO 1 I=1,NX
         X(I)=DFLOAT(I-1)/DFLOAT(NX-1)
1     CONTINUE
C...
C...  INITIAL CONDITIONS
      DO 2 I=1,NX
         U1(I)=   DSIN(C3*X(I))
         U2(I)=C1*DSIN(C3*X(I))
2     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV1
      IP=0
      RETURN
      END
      SUBROUTINE INIT2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/       T,   NSTOP,   NORUN
     1      /Y/  U1(NX),  U2(NX)
     2      /F/ U1T(NX), U2T(NX)
     3      /S/ U1X(NX),U1XX(NX),   X(NX)
     4      /C/       A,       B,       C,       D,
     5               C1,      C2,      PI
     6      /I/      IP
C...
C...  PROBLEM PARAMETERS
      A= 1.0D0
      B= 1.0D0
      C= 1.0D0
      D=-1.0D0
C...
C...  WEAK NONLINEAR DAMPING
      IF(NORUN.EQ.3)THEN
         C2=1.0D0
C...
C...  STRONG NONLINEAR DAMPING
      ELSE IF(NORUN.EQ.4)THEN
         C2=100.0D0
      END IF
C...
C...  PRECOMPUTE SOME CONSTANTS FOR USE IN THE INITIAL CONDITIONS
      C3=4.0D0*DATAN(1.0D0)
      C1=-B/(2.0D0*A)
C...
C...  SPATIAL GRID
      DO 1 I=1,NX
         X(I)=DFLOAT(I-1)/DFLOAT(NX-1)
1     CONTINUE
C...
C...  INITIAL CONDITIONS
      DO 2 I=1,NX
         U1(I)=   DSIN(C3*X(I))
         U2(I)=C1*DSIN(C3*X(I))
2     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV2
      IP=0
      RETURN
      END
      SUBROUTINE DERV1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/       T,   NSTOP,   NORUN
     1      /Y/  U1(NX),  U2(NX)
     2      /F/ U1T(NX), U2T(NX)
     3      /S/ U1X(NX),U1XX(NX),   X(NX)
     4      /C/       A,       B,       C,       D,
     5               C1,      C2,      C3
     6      /I/      IP
C...
C...  BOUNDARY CONDITION AT X = 0
      NL=1
      U1(1)=0.0D0
C...
C...  BOUNDARY CONDITION AT X = 1
      NU=1
      U1(NX)=0.0D0
C...
C...  UXX
      CALL DSS044(0.0D0,1.0D0,NX,U1,U1X,U1XX,NL,NU)
C...
C...  PDES
      DO 1 I=1,NX
      IF((I.EQ.1).OR.(I.EQ.NX))THEN
         U1T(I)=0.0D0
      ELSE
         U1T(I)=U2(I)
      END IF
         U2T(I)=(1.0D0/A)*(C*U1XX(I)-B*U1T(I)+D*U1(I))
1     CONTINUE
      RETURN
      END
      SUBROUTINE DERV2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/       T,   NSTOP,   NORUN
     1      /Y/  U1(NX),  U2(NX)
     2      /F/ U1T(NX), U2T(NX)
     3      /S/ U1X(NX),U1XX(NX),   X(NX)
     4      /C/       A,       B,       C,       D,
     5               C1,      C2,      PI
     6      /I/      IP
C...
C...  BOUNDARY CONDITION AT X = 0
      NL=1
      U1(1)=0.0D0
C...
C...  BOUNDARY CONDITION AT X = 1
      NU=1
      U1(NX)=0.0D0
C...
C...  UXX
      CALL DSS044(0.0D0,1.0D0,NX,U1,U1X,U1XX,NL,NU)
C...
C...  PDES
      DO 1 I=1,NX
      IF((I.EQ.1).OR.(I.EQ.NX))THEN
         U1T(I)=0.0D0
      ELSE
         U1T(I)=U2(I)
      END IF
         U2T(I)=(1.0D0/A)*(C*U1XX(I)-B*U1T(I)+D*U1(I))
     1         -C2*DABS(U1T(I))*U1T(I)
1     CONTINUE
      RETURN
      END
      DOUBLE PRECISION FUNCTION UANAL(X,T)
C...
C...  FUNCTION UANAL COMPUTES THE ANALYTICAL SOLUTION
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/C/       A,       B,       C,       D,
     5               C1,      C2,      C3
      UANAL=DEXP(C1*T)*DCOS(C2*T)*DSIN(C3*X)
      RETURN
      END
      SUBROUTINE PRINT1(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/       T,   NSTOP,   NORUN
     1      /Y/  U1(NX),  U2(NX)
     2      /F/ U1T(NX), U2T(NX)
     3      /S/ U1X(NX),U1XX(NX),   X(NX)
     4      /C/       A,       B,       C,       D,
     5               C1,      C2,      C3
     6      /I/      IP
      DIMENSION  UA(NX)
C...
C...  INITIALIZE A COUNTER FOR THE PRINTING AND PLOTTING
      DATA IP/0/
      IP=IP+1
      IF(IP.EQ.1)THEN
C...
C...  MAP THE ODE JACOBIAN MATRIX AT THE BEGINNING OF THE SOLUTION
      CALL MAP
      END IF
C...
C...  MONITOR THE CALCULATIONS
      WRITE(*,*)NORUN,T
C...
C...  COMPUTE THE ANALYTICAL SOLUTION AT X = 0, 0.2, ..., 1
      DO 1 I=1,NX,4
         UA(I)=UANAL(X(I),T)
1     CONTINUE
C...
C...  PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS AT X = 0, 0.2,
C...  ..., 1
      IF(T.EQ.0.)WRITE(NO,4)
4     FORMAT(10X,'T',8X,'X=0',3X,'   X=0.2',3X,'   X=0.4',
     1                        3X,'   X=0.6',3X,'   X=0.8',
     2                        7X,' X=1')
      WRITE(NO,2)T,(U1(I),I=1,NX,4)
2     FORMAT(F11.2,/,'  NUMERICAL',6F11.5)
      WRITE(NO,3)(UA(I),I=1,NX,4)
3     FORMAT(' ANALYTICAL',6F11.5,/)
C...
C...  MAP THE ODE JACOBIAN MATRIX AT THE END OF THE SOLUTION
      IF(IP.LT.11)RETURN
      CALL MAP
C...
C...  RESET THE INTEGER COUNTER FOR THE NEXT RUN
      IP=0
      RETURN
      END
      SUBROUTINE PRINT2(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/       T,   NSTOP,   NORUN
     1      /Y/  U1(NX),  U2(NX)
     2      /F/ U1T(NX), U2T(NX)
     3      /S/ U1X(NX),U1XX(NX),   X(NX)
     4      /C/       A,       B,       C,       D,
     5               C1,      C2,      PI
     6      /I/      IP
C...
C...  INITIALIZE A COUNTER FOR THE PRINTING AND PLOTTING
      DATA IP/0/
      IP=IP+1
      IF(IP.EQ.1)THEN
C...
C...  MAP THE ODE JACOBIAN MATRIX AT THE BEGINNING OF THE SOLUTION
      CALL MAP
      END IF
C...
C...  MONITOR THE CALCULATIONS
      WRITE(*,*)NORUN,T
C...
C...  PRINT THE NUMERICAL SOLUTION AT X = 0, 0.2,..., 1
      IF(T.EQ.0.)WRITE(NO,4)
4     FORMAT(10X,'T',8X,'X=0',3X,'   X=0.2',3X,'   X=0.4',
     1                        3X,'   X=0.6',3X,'   X=0.8',
     2                        7X,' X=1')
      WRITE(NO,2)T,(U1(I),I=1,NX,4)
2     FORMAT(F11.2,/,'  NUMERICAL',6F11.5,/)
C...
C...  MAP THE ODE JACOBIAN MATRIX AT THE END OF THE SOLUTION
      IF(IP.LT.11)RETURN
      CALL MAP
C...
C...  RESET THE INTEGER COUNTER FOR THE NEXT RUN
      IP=0
      RETURN
      END
      SUBROUTINE MAP
      PARAMETER (N=42)
      COMMON/T/      T
     1      /Y/   Y(N)
     2      /F/   F(N)
C...
C...  SUBROUTINE MAP CALLS SUBROUTINE JMAP TO MAP THE JACOBIAN MATRIX 
C...  OF AN NTH-ORDER ODE SYSTEM.
C...
C...  DECLARE SELECTED VARIABLES DOUBLE PRECISION
      DOUBLE PRECISION T, Y, F, YOLD(N), FOLD(N)
C...
C...  DEFINE SINGLE PRECISION ARRAY REQUIRED BY SUBROUTINE JMAP (A)
      REAL A(N,N)
C...
C...  MAP THE JACOBIAN MATRIX OF THE ODE SYSTEM DEFINED IN SUBROUTINE
C...  DERV
      CALL JMAP(N,A,Y,YOLD,F,FOLD)
      RETURN
      END

