      SUBROUTINE INITAL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/         T,     NSTOP,     NORUN

C...  V(0,T) = 1 - EXP(-C2*T)
      IF((NORUN-1)*(NORUN-2)*(NORUN-3).EQ.0)CALL INIT1      
C...
C...  V(0,T) = EXP(-C2*T)
      IF((NORUN-4)*(NORUN-5)*(NORUN-6).EQ.0)CALL INIT2      
C...
C...  V(0,T) = 1 (UNIT STEP)
      IF((NORUN-7)*(NORUN-8)*(NORUN-9).EQ.0)CALL INIT3      
C...
C...  V(0,T) = T (UNIT RAMP)
      IF((NORUN-10)*(NORUN-11)*(NORUN-12).EQ.0)CALL INIT4      
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/         T,     NSTOP,     NORUN
C...
C...  V(0,T) = 1 - EXP(-C2*T)
      IF((NORUN-1)*(NORUN-2)*(NORUN-3).EQ.0)CALL DERV1
C...
C...  V(0,T) = EXP(-C2*T)
      IF((NORUN-4)*(NORUN-5)*(NORUN-6).EQ.0)CALL DERV2     
C...
C...  V(0,T) = 1 (UNIT STEP)
      IF((NORUN-7)*(NORUN-8)*(NORUN-9).EQ.0)CALL DERV3     
C...
C...  V(0,T) = T (UNIT RAMP)
      IF((NORUN-10)*(NORUN-11)*(NORUN-12).EQ.0)CALL DERV4      
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/T/         T,     NSTOP,     NORUN
C...
C...  V(0,T) = 1 - EXP(-C2*T)
      IF((NORUN-1)*(NORUN-2)*(NORUN-3).EQ.0)CALL PRINT1(NI,NO)
C...
C...  V(0,T) = EXP(-C2*T)
      IF((NORUN-4)*(NORUN-5)*(NORUN-6).EQ.0)CALL PRINT2(NI,NO)
C...
C...  V(0,T) = 1 (UNIT STEP)
      IF((NORUN-7)*(NORUN-8)*(NORUN-9).EQ.0)CALL PRINT3(NI,NO)
C...
C...  V(0,T) = T (UNIT RAMP)
      IF((NORUN-10)*(NORUN-11)*(NORUN-12).EQ.0)CALL PRINT4(NI,NO)      
      RETURN
      END
      SUBROUTINE INIT1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/    U1(NX),    U2(NX),       U1A
     +      /F/   U1T(NX),   U2T(NX),      U1AT
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,       ARG
C...
C...  SET THE PROBLEM PARAMETERS
         C1=0.560D+00
         C2=0.100D+00
          V=2.031D+00
          X=1.0D+00
C...
C...  PRECOMPUTE SOME QUANTITIES USED IN THE ANALYTICAL SOLUTION
      EXP1=DEXP(-C1/V*X)
C...
C...  INITIAL CONDITIONS
      DO 1 I=1,NX
         U1(I)=0.0D+00
         U2(I)=0.0D+00
1     CONTINUE
      U1A=0.0D+00
C...
C...  INITIAL DERIVATIVES
      CALL DERV1
      RETURN
      END
      SUBROUTINE INIT2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/    U1(NX),    U2(NX),       U1A
     +      /F/   U1T(NX),   U2T(NX),      U1AT
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,       ARG
C...
C...  SET THE PROBLEM PARAMETERS
         C1=0.560D+00
         C2=0.100D+00
          V=2.031D+00
          X=1.0D+00
C...
C...  PRECOMPUTE SOME QUANTITIES USED IN THE ANALYTICAL SOLUTION
      EXP1=DEXP(-C1/V*X)
C...
C...  INITIAL CONDITIONS
      DO 1 I=1,NX
         U1(I)=0.0D+00
         U2(I)=0.0D+00
1     CONTINUE
      U1A=0.0D+00
C...
C...  INITIAL DERIVATIVES
      CALL DERV2
      RETURN
      END
      SUBROUTINE INIT3
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/    U1(NX),    U2(NX),       U1A
     +      /F/   U1T(NX),   U2T(NX),      U1AT
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,       ARG
C...
C...  SET THE PROBLEM PARAMETERS
         C1=0.560D+00
         C2=0.100D+00
          V=2.031D+00
          X=1.0D+00
C...
C...  PRECOMPUTE SOME QUANTITIES USED IN THE ANALYTICAL SOLUTION
      EXP1=DEXP(-C1/V*X)
C...
C...  INITIAL CONDITIONS
      DO 1 I=1,NX
         U1(I)=0.0D+00
         U2(I)=0.0D+00
1     CONTINUE
      U1A=0.0D+00
C...
C...  INITIAL DERIVATIVES
      CALL DERV3
      RETURN
      END
      SUBROUTINE INIT4
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21,NT=3)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/    U1(NX),    U2(NX),   U1A(NT)
     +      /F/   U1T(NX),   U2T(NX),  U1AT(NT)
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,       ARG
C...
C...  SET THE PROBLEM PARAMETERS
         C1=0.560D+00
         C2=0.100D+00
          V=2.031D+00
          X=1.0D+00
C...
C...  PRECOMPUTE SOME QUANTITIES USED IN THE ANALYTICAL SOLUTION
      EXP1=DEXP(-C1/V*X)
C...
C...  INITIAL CONDITIONS
      DO 1 I=1,NX
         U1(I)=0.0D+00
         U2(I)=0.0D+00
1     CONTINUE
      DO 2 I=1,NT
         U1A(I)=0.0D+00
2     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV4
      RETURN
      END
      SUBROUTINE DERV1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/    U1(NX),    U2(NX),       U1A
     +      /F/   U1T(NX),   U2T(NX),      U1AT
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,   ARG
C...
C...  BOUNDARY CONDITION
       U1(1)=1.0D+00-DEXP(-C2*T)
      U1T(1)=0.0D+00
C...
C...  U1X
C...
C...     TWO POINT UPWIND
         IF(NORUN.EQ.1)CALL DSS012(0.D+0,X,NX,U1,U1X,1.0D+00)
C...
C...     FOUR POINT BIASED UPWIND
         IF(NORUN.EQ.2)CALL DSS018(0.D+0,X,NX,U1,U1X,1.0D+00)
C...
C...     FIVE POINT BIASED UPWIND
         IF(NORUN.EQ.3)CALL DSS020(0.D+0,X,NX,U1,U1X,1.0D+00)
C...
C...  PDES
      DO 1 I=1,NX
         IF(I.NE.1)THEN
            U1T(I)=C1*(U2(I)-U1(I))-V*U1X(I)
         END IF
            U2T(I)=C2*(U1(I)-U2(I))
1     CONTINUE
C...
C...  TEST FOR THE ELAPSED TIME
C...  T - X/V LT 0
      IF((T-X/V).LT.0.0D+00)THEN
       U1AT=0.0D+00
C...
C...  T - X/V GE 0
      ELSE
         ARG=2.0D+00*DSQRT((C1*C2*X/V)*(T-X/V))
        U1AT=EXP1*C2*DEXP(-C2*(T-X/V))*BESSI0(ARG)
      END IF
      RETURN
      END
      SUBROUTINE DERV2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  DECLARE SELECTED VARIABLES AS DOUBLE PRECISION
      DOUBLE PRECISION  K,     V,     C,     X,     U,     
     +               EXP1,   ARG,  U1AT,   U1A,     T,
     +                 U1,    U2,   U1T,   U2T,   U1X,
     +             BESSI0
      PARAMETER (NX=21)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/    U1(NX),    U2(NX),       U1A
     +      /F/   U1T(NX),   U2T(NX),      U1AT
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,   ARG
C...
C...  BOUNDARY CONDITION
       U1(1)=DEXP(-C2*T)
      U1T(1)=0.0D+00
C...
C...  U1X
C...
C...     TWO POINT UPWIND
         IF(NORUN.EQ.4)CALL DSS012(0.0D+00,X,NX,U1,U1X,1.0D+00)
C...
C...     FOUR POINT BIASED UPWIND
         IF(NORUN.EQ.5)CALL DSS018(0.0D+00,X,NX,U1,U1X,1.0D+00)
C...
C...     FIVE POINT BIASED UPWIND
         IF(NORUN.EQ.6)CALL DSS020(0.0D+00,X,NX,U1,U1X,1.0D+00)
C...
C...  PDES
      DO 1 I=1,NX
         IF(I.NE.1)THEN
            U1T(I)=C1*(U2(I)-U1(I))-V*U1X(I)
         END IF
            U2T(I)=C2*(U1(I)-U2(I))
1     CONTINUE
C...
C...  TEST FOR THE ELAPSED TIME
C...  T - X/V LT 0
      IF((T-X/V).LT.0.0D+00)THEN
       U1AT=0.0D+00
C...
C...  T - X/V GE 0
      ELSE
        ARG=2.0D+00*DSQRT((C1*C2*X/V)*(T-X/V))
       U1AT=EXP1*DEXP(-C2*(T-X/V))*BESSI0(ARG)
      END IF
      RETURN
      END
      SUBROUTINE DERV3
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/    U1(NX),    U2(NX),       U1A
     +      /F/   U1T(NX),   U2T(NX),      U1AT
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,   ARG
C...
C...  BOUNDARY CONDITION
       U1(1)=1.0D+00
      U1T(1)=0.0D+00
C...
C...  U1X
C...
C...     TWO POINT UPWIND
         IF(NORUN.EQ.7)CALL DSS012(0.0D+00,X,NX,U1,U1X,1.0D+00)
C...
C...     FOUR POINT BIASED UPWIND
         IF(NORUN.EQ.8)CALL DSS018(0.0D+00,X,NX,U1,U1X,1.0D+00)
C...
C...     FIVE POINT BIASED UPWIND
         IF(NORUN.EQ.9)CALL DSS020(0.0D+00,X,NX,U1,U1X,1.0D+00)
C...
C...  PDES
      DO 1 I=1,NX
         IF(I.NE.1)THEN
            U1T(I)=C1*(U2(I)-U1(I))-V*U1X(I)
         END IF
            U2T(I)=C2*(U1(I)-U2(I))
1     CONTINUE
C...
C...  TEST FOR THE ELAPSED TIME
C...  T - X/V LT 0
      IF((T-X/V).LT.0.0D+00)THEN
       U1AT=0.0D+00
C...
C...  T - X/V GE 0
      ELSE
        ARG=2.0D+00*DSQRT((C1*C2*X/V)*(T-X/V))
       U1AT=EXP1*C2*DEXP(-C2*(T-X/V))*BESSI0(ARG)
      END IF
      RETURN
      END      
      SUBROUTINE DERV4
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21,NT=3)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/    U1(NX),    U2(NX),   U1A(NT)
     +      /F/   U1T(NX),   U2T(NX),  U1AT(NT)
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,   ARG
C...
C...  BOUNDARY CONDITION
      IF(T.LE.1.0D+00)THEN
         U1(1)=T
        U1T(1)=1.0D+00
      ELSE IF(T.GT.1.0D+00)THEN
         U1(1)=1.0D+00
        U1T(1)=0.0D+00
      END IF
C...
C...  U1X
C...
C...     TWO POINT UPWIND
         IF(NORUN.EQ.10)CALL DSS012(0.0D+00,X,NX,U1,U1X,1.0D+00)
C...
C...     FOUR POINT BIASED UPWIND
         IF(NORUN.EQ.11)CALL DSS018(0.0D+00,X,NX,U1,U1X,1.0D+00)
C...
C...     FIVE POINT BIASED UPWIND
         IF(NORUN.EQ.12)CALL DSS020(0.0D+00,X,NX,U1,U1X,1.0D+00)
C...
C...  PDES
      DO 1 I=1,NX
         IF(I.NE.1)THEN
            U1T(I)=C1*(U2(I)-U1(I))-V*U1X(I)
         END IF
            U2T(I)=C2*(U1(I)-U2(I))
1     CONTINUE
C...
C...  TEST FOR THE ELAPSED TIME FOR RESPONSE TO STEP AT T = 0
C...  T - X/V LT 0
      IF((T-X/V).LT.0.0D+00)THEN
       U1AT(1)=0.0D+00
C...
C...  T - X/V GE 0
      ELSE
         ARG=2.0D+00*DSQRT((C1*C2*X/V)*(T-X/V))
       U1AT(1)=EXP1*C2*DEXP(-C2*(T-X/V))*BESSI0(ARG)
      END IF
C...
C...  TEST FOR THE ELAPSED TIME FOR RESPONSE TO STEP AT T = 1
C...  T - X/V - 1 LT 0
      IF(((T-X/V)-1.0D+00).LT.0.0D+00)THEN
       U1AT(2)=0.0D+00
C...
C...  T - X/V - 1 GE 0
      ELSE
          ARG=2.0D+00*DSQRT((C1*C2*X/V)*(T-X/V-1.0D+00))
      U1AT(2)=-EXP1*C2*DEXP(-C2*(T-X/V-1.0D+00))*BESSI0(ARG)
      END IF
C...
C...  SUPERIMPOSE TWO STEP RESPONSES
      U1AT(3)=U1AT(1)/C2+U1A(1)+U1AT(2)/C2+U1A(2)
      RETURN
      END
      SUBROUTINE PRINT1(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/    U1(NX),    U2(NX),       U1A
     +      /F/   U1T(NX),   U2T(NX),      U1AT
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,   ARG
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
C...
C...  PRINT A HEADING
      WRITE(NO,1)
1     FORMAT(//,14X,'T',8X,'U1(X,T)',3X,'U1(X,T) ANAL')
      END IF
C...
C...  PRINT THE SOLUTION 
      WRITE(NO,2)T,U1(NX),U1A
2     FORMAT(F15.2,2F15.4)
C...
C...  WRITE NORUN, IP AND T ON THE SCREEN TO MONITOR THE PROGRESS OF
C...  THE CALCULATIONS
      WRITE(*,*)' NORUN = ',NORUN,' IP = ',IP,' T = ',T
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
     +      /Y/    U1(NX),    U2(NX),       U1A
     +      /F/   U1T(NX),   U2T(NX),      U1AT
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,       ARG
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
C...  CALL MAP
C...
C...  PRINT A HEADING
      WRITE(NO,1)
1     FORMAT(//,14X,'T',8X,'U1(X,T)',3X,'U1(X,T) ANAL')
      END IF
C...
C...  PRINT THE SOLUTION
      WRITE(NO,2)T,U1(NX),U1AT
2     FORMAT(F15.2,2F15.4)
C...
C...  WRITE NORUN, IP AND T ON THE SCREEN TO MONITOR THE PROGRESS OF
C...  THE CALCULATIONS
      WRITE(*,*)' NORUN = ',NORUN,' IP = ',IP,' T = ',T
C...
C...  MAP THE ODE JACOBIAN MATRIX AT THE END OF THE SOLUTION
      IF(IP.LT.11)RETURN
C...  CALL MAP
C...
C...  RESET THE INTEGER COUNTER FOR THE NEXT RUN
      IP=0
      RETURN
      END
      SUBROUTINE PRINT3(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=21)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/    U1(NX),    U2(NX),       U1A
     +      /F/   U1T(NX),   U2T(NX),      U1AT
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,       ARG
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
C...  CALL MAP
C...
C...  PRINT A HEADING
      WRITE(NO,1)
1     FORMAT(//,14X,'T',8X,'U1(X,T)',3X,'U1(X,T) ANAL')
      END IF
C...
C...  PRINT THE SOLUTION
      WRITE(NO,2)T,U1(NX),U1AT/C2+U1A
2     FORMAT(F15.2,2F15.4)
C...
C...  WRITE NORUN, IP AND T ON THE SCREEN TO MONITOR THE PROGRESS OF
C...  THE CALCULATIONS
      WRITE(*,*)' NORUN = ',NORUN,' IP = ',IP,' T = ',T
C...
C...  MAP THE ODE JACOBIAN MATRIX AT THE END OF THE SOLUTION
      IF(IP.LT.11)RETURN
C...  CALL MAP
C...
C...  RESET THE INTEGER COUNTER FOR THE NEXT RUN
      IP=0
      RETURN
      END      
      SUBROUTINE PRINT4(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  DECLARE SELECTED VARIABLES AS DOUBLE PRECISION
      PARAMETER (NX=21,NT=3)
      COMMON/T/         T,     NSTOP,     NORUN
     +      /Y/    U1(NX),    U2(NX),   U1A(NT)
     +      /F/   U1T(NX),   U2T(NX),  U1AT(NT)
     +      /S/   U1X(NX)
     +      /C/        C1,        C2,         V,         X,
     +               EXP1,       ARG
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
C...  CALL MAP
C...
C...  PRINT A HEADING
      WRITE(NO,1)
1     FORMAT(//,14X,'T',8X,'U1(X,T)',3X,'U1(X,T) ANAL')
      END IF
C...
C...  PRINT THE SOLUTION
      WRITE(NO,2)T,U1(NX),U1A(3)
2     FORMAT(F15.2,2F15.4)
C...
C...  WRITE NORUN, IP AND T ON THE SCREEN TO MONITOR THE PROGRESS OF
C...  THE CALCULATIONS
      WRITE(*,*)' NORUN = ',NORUN,' IP = ',IP,' T = ',T
C...
C...  MAP THE ODE JACOBIAN MATRIX AT THE END OF THE SOLUTION
      IF(IP.LT.11)RETURN
C...  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
      DOUBLE PRECISION FUNCTION BESSI0(X)
C...  FUNCTION BESSI0(X)
C...
C...  FUNCTION BESSI0 IS A MINOR REVISION BY W. E. SCHIESSER OF THE
C...  ORIGINAL NUMERICAL RECIPES FUNCTION TO BRING IT INTO CLOSER
C...  CONFORMITY WITH STANDARD FORTRAN DOUBLE PRECISION PROGRAMMING.
C...  THE ORIGINAL STATEMENTS THAT HAVE BEEN MODIFIED ARE COMMENTED
C...  OUT
      DOUBLE PRECISION     X,    AX,
     *       Y,    P1,    P2,    P3,    P4,   P5,    P6,    P7,
     *      Q1,    Q2,    Q3,    Q4,    Q5,   Q6,    Q7,    Q8,
     *      Q9
C...  REAL*8 Y,    P1,    P2,    P3,    P4,   P5,    P6,    P7,
C... *      Q1,    Q2,    Q3,    Q4,    Q5,   Q6,    Q7,    Q8,
C... *      Q9
      DATA         P1,           P2,          P3,           P4,
     *             P5,           P6,          P7/
     *          1.0D0,  3.5156229D0, 3.0899424D0,  1.2067492D0,
     *    0.2659732D0,  0.360768D-1,  0.45813D-2/
      DATA         Q1,           Q2,          Q3,           Q4,
     *             Q5,           Q6,          Q7,           Q8,
     *             Q9/
     *   0.39894228D0, 0.1328592D-1, 0.225319D-2, -0.157565D-2,
     *    0.916281D-2,-0.2057706D-1,0.2635537D-1,-0.1647633D-1,
     *    0.392377D-2/
      IF (DABS(X).LT.3.75D0) THEN
C...  IF ( ABS(X).LT.3.75  ) THEN
        Y=(X/3.75D0)**2
C...    Y=(X/3.75  )**2
        BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
      ELSE
        AX=DABS(X)
C...    AX= ABS(X)
        Y=3.75D0/AX
C...    Y=3.75  /AX
        BESSI0=(DEXP(AX)/DSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4
     *      +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
C...    BESSI0=( EXP(AX)/ SQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4
C... *      +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
      ENDIF
      RETURN
      END
