C***********************************************************************

      SUBROUTINE INITIAL(DT,DX,Nx,Q,Delta,acon,Ndir)
      
      IMPLICIT NONE

      INTEGER j, Nx, Nnon0, Nfront, Nback, Ndir
      REAL*8  DT, DX, Q(-1:1,-Nx:Nx), Delta, Qdum, front, back, test, 
     /        acon

C Set the initial conditions of the Solution

      DO j = -Nx, Nx
         Q(1,j) = 0.0D0
      ENDDO 
      
      Nnon0 = INT(Nx*0.2D0)
      DO j = -Nnon0, Nnon0
         Q(1,j) = 0.04D0*EXP( -(DX*j)**2 / Delta**2 )
      ENDDO 

C Now determine your backwards step using the analytical solution.

 400  DO j = -Nx, Nx
         Q(0,j) = 0.0D0
      ENDDO

      back = ( -0.2D0 - Ndir*DT )*EXP( acon*DT )
      Nback =  INT(back/DX)
      test = ABS( NINT(back/DX) - back/DX )
      IF ( test .LE. 1.0D-7 ) Nback = NINT(back/DX)

      front = ( 0.2D0 - Ndir*DT )*EXP( acon*DT )
      Nfront = INT(front/DX)
      test = ABS( NINT(front/DX) - front/DX )
      IF ( test .LE. 1.0D-7 ) Nfront = NINT(front/DX)

      DO j = Nback, Nfront
         Q(0,j) = 0.04D0*EXP( -( DX*j*EXP( -acon*DT ) + Ndir*DT )**2 
     /                        / Delta**2 )
      ENDDO 

      RETURN 
      END

C***********************************************************************

      SUBROUTINE MATRIX(DT,DX,itime,thetax,thetat,thetaf,acon,Nbound,
     /                  Nuse,Mat,Work,Ord)

      IMPLICIT NONE

      INTEGER itime, j, k, Nuse, Nbound, info
      REAL*8  DT, DX, rho, thetax, thetat, thetaf, acon, beta, gamma,
     /        Mat(Nuse,Nuse), Work(Nuse), rcond, thxa, thxb, gu11,
     /        sqrtgu, test
      INTEGER Ord(Nuse)

C Set the initial conditions of the matrix.

      DO j = 1, Nuse
         Work(j) = 0.0D0
         DO k = 1, Nuse
            Mat(j,k) = 0.0D0
         ENDDO
      ENDDO

      rho = DT / DX
      gu11 = EXP( -2.0D0*acon*(itime*DT) )
      sqrtgu = EXP( -acon*( (itime+1)*DT ) )
      IF ( acon .LT. 0.0D0 ) THEN
         test = ( -1.0D0 + (itime+1)*DT )*sqrtgu
      ELSE
         test = ( 1.0D0 + (itime+1)*DT )*sqrtgu
      ENDIF

      IF ( Nbound .EQ. 1 ) THEN
         IF ( test .LT. 1.0D0 .AND. acon .GE. 0.0D0 ) THEN
            Mat(1,1) = 1.0D0 
         ELSE
            Mat(1,1) = 3.0D0*( rho*( sqrtgu - acon ) + 1.0D0 )
            Mat(1,2) = -4.0D0*rho*( sqrtgu - acon )
            Mat(1,3) = rho*( sqrtgu - acon )
         ENDIF
         DO j = 2, Nuse-1
            beta = ( ( j - 1.0D0 )*DX - 1.0D0 )*acon
            gamma = ( ( j - 1.0D0 )*DX - 1.0D0 )*acon**2
            thxa = ( gu11 - beta**2 )*thetax*rho**2
            thxb = rho*( beta - 0.5D0*gamma*DT*thetaf )

            Mat(j,j-1) = 0.5D0*( thetat - thxa + thxb )
            Mat(j,j) = 1.0D0 - thetat + thxa
            Mat(j,j+1) = 0.5D0*( thetat - thxa - thxb )
         ENDDO
         IF ( test .LT. 1.0D0 .AND. acon .GE. 0.0D0 ) THEN
            Mat(Nuse,Nuse) = 1.0D0
         ELSE
            Mat(Nuse,Nuse-2) = rho*( sqrtgu - acon )
            Mat(Nuse,Nuse-1) = -4.0D0*rho*( sqrtgu - acon )
            Mat(Nuse,Nuse) = 3.0D0*( rho*( sqrtgu - acon ) + 1.0D0 )
         ENDIF
      ELSE
         PRINT *, 'Invalid value of Nbound' 
         STOP
      ENDIF

C     CALL DGEFA(Mat,Nuse,Nuse,Ord,INFO)
C     IF ( INFO .NE. 0 ) THEN
C        PRINT *, 'The Inversion has failed.'
C        STOP
C     ENDIF
      CALL DGECO(Mat,Nuse,Nuse,Ord,rcond,work)

      RETURN
      END

C***********************************************************************

      SUBROUTINE EVOL(DT,DX,itime,thetax,thetat,thetaf,acon,Nbound,Nx,Q,
     /                Nuse,Mat,Vec,Ord)

      IMPLICIT NONE

      INTEGER j, k, itime, Nx, Nuse, Nbound
      REAL*8  DT, DX, thetax, thetat, thetaf, thxa, thxb, thxc, thxd,
     /        Q(-1:1,-Nx:Nx), Mat(Nuse,Nuse), Vec(Nuse), acon, beta, 
     /        gamma, rho, gu11, test
      INTEGER Ord(Nuse)

C Evolve to the next timestep.

      DO j = -Nx, Nx
         Q(-1,j) = Q(0,j)
         Q(0,j) = Q(1,j)
      ENDDO

      rho = ( DT / DX )
      gu11 = EXP( -2.0*acon*(itime*DT) )
      IF ( acon .GE. 0.0D0 )
     /    test = ( 1.0D0 + (itime+1)*DT )*EXP( -acon*( (itime+1)*DT ) )

      DO j = -Nx+1, Nx-1
         beta = j*DX*acon
         gamma = j*DX*acon**2
         thxa = thetax*( gu11 - beta**2 )*rho**2
         thxb = ( 1.0D0 - thetax )*( gu11 - beta**2 )*rho**2
         thxc = 0.5D0*gamma*rho*DT*( 1.0D0 - thetaf )
         thxd = 0.5D0*rho*( beta + 0.5D0*thetaf*gamma*DT )
         k = j + Nx + 1
         Vec(k) = ( thxb + thetat )*( Q(0,j+1) + Q(0,j-1) ) 
     /            - thxc*( Q(0,j+1) - Q(0,j-1) ) 
     /            + 2.0D0*( 1.0D0 - thetat - thxb )*Q(0,j)
     /            - ( 1.0D0 - thetat + thxa )*Q(-1,j)
     /            + 0.5D0*( thxa - thetat )*( Q(-1,j+1) + Q(-1,j-1) )
     /            - thxd*( Q(-1,j+1) - Q(-1,j-1) )
      ENDDO
      IF ( test .LT. 1.0D0 .AND. acon .GE. 0.0D0 ) THEN
         Vec(1) = 0.0D0
         Vec(Nuse) = 0.0D0
      ELSE
         Vec(1) = 4.0D0*Q(0,-Nx) - Q(-1,-Nx) 
         Vec(Nuse) = 4.0D0*Q(0,Nx) - Q(-1,Nx)  
      ENDIF

      CALL DGESL(Mat,Nuse,Nuse,Ord,Vec,0)

      DO j = -Nx, Nx
         k = j + Nx + 1
         Q(1,j) = Vec(k)
      ENDDO

      RETURN

      END