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

      SUBROUTINE INITIAL(DT,DX,Nx,Q,Delta,beta,Cloak)
      
      IMPLICIT NONE

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

C Set the initial conditions of the Solution

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

C Now determine your backwards step using the analytical solution.

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

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

      front = 0.2D0 - Sdt
      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 + Sdt)**2 / Delta**2 )
      ENDDO 

      RETURN 
      END

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

      SUBROUTINE MATRIX(DT,DX,thetax,thetat,beta,Nbound,Nuse,Mat,Work,
     /                  Ord)

      IMPLICIT NONE

      INTEGER j, k, Nuse, Nbound, info
      REAL*8  DT, DX, rho, thetax, thetat, beta, Mat(Nuse,Nuse), 
     /        Work(Nuse), rcond
      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

      IF ( Nbound .EQ. -1 ) THEN
         DO j = 1, Nuse
            IF ( j .NE. 1 ) Mat(j,j-1) = -rho**2
            Mat(j,j) = 1.0D0 + 2.0D0*rho**2
            IF ( j .NE. Nuse ) Mat(j,j+1) = -rho**2
         ENDDO
      ELSEIF ( Nbound .EQ. 0 ) THEN
         Mat(1,Nuse) = -rho**2
         Mat(1,1) = 1.0D0 + 2.0D0*rho**2
         Mat(1,2) = -rho**2
         DO j = 2, Nuse
            Mat(j,j-1) = -rho**2
            Mat(j,j) = 1.0D0 + 2.0D0*rho**2
            IF ( j .EQ. Nuse ) THEN
               Mat(j,1) = -rho**2
            ELSE
              Mat(j,j+1) = -rho**2
            ENDIF
         ENDDO
      ELSEIF ( Nbound .EQ. 1 ) THEN
         Mat(1,1) = 3.0D0*( rho + 1.0D0 )
         Mat(1,2) = -4.0D0*rho
         Mat(1,3) =  rho
         DO j = 2, Nuse-1
            Mat(j,j-1) = -rho**2
            Mat(j,j) = 1.0D0 + 2.0D0*rho**2
            Mat(j,j+1) = -rho**2
         ENDDO
         Mat(Nuse,Nuse-2) = rho
         Mat(Nuse,Nuse-1) = -4.0D0*rho
         Mat(Nuse,Nuse) = 3.0D0*( rho + 1.0D0 )
      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)
      PRINT 15, 1.0D0 / rcond
 15   FORMAT(' The condition number = ',G15.7)

      RETURN
      END



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

      SUBROUTINE EVOL(DT,DX,thetax,thetat,beta,Nbound,Nx,Q,Nuse,Mat,
     /                Vec,Ord,Cloak,Ncloak)

      IMPLICIT NONE

      INTEGER j, k, Nx, Nuse, Nbound, Ncloak
      REAL*8  DT, DX, thetax, thetat, beta, rho, Q(-1:1,-Nx:Nx), 
     /        Mat(Nuse,Nuse), Vec(Nuse), Cloak(-Nx:Nx)
      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 )

      IF ( Nbound .EQ. -1 ) THEN
         DO j = -Nx+1, Nx-1
            k = j + Nx
            Vec(k) = 2.0D0*Q(0,j) - Q(-1,j)
         ENDDO
      ELSEIF ( Nbound .EQ. 0 ) THEN
         DO j = -Nx, Nx-1
            k = j + Nx + 1
            Vec(k) = 2.0D0*Q(0,j) - Q(-1,j)
         ENDDO
      ELSE
         Vec(1) = 4.0D0*Q(0,-Nx) - Q(-1,-Nx)
         DO j = -Nx+1, Nx-1
            k = j + Nx + 1
            Vec(k) = 2.0D0*Q(0,j) - Q(-1,j)
         ENDDO
         Vec(Nuse) = 4.0D0*Q(0,Nx) - Q(-1,Nx)   
      ENDIF

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

      IF ( Nbound .EQ. -1 ) THEN
         Q(1,-Nx) = 0.0D0      
         DO j = -Nx+1, Nx-1
            k = j + Nx
            Q(1,j) = Vec(k)
         ENDDO
         Q(1,Nx) = 0.0D0
      ELSE
         DO j = -Nx, Nx-1
            k = j + Nx + 1
            Q(1,j) = Vec(k)
         ENDDO
         IF ( Nbound .EQ. 0 ) THEN
            Q(1,Nx) = Q(1,-Nx)
         ELSE 
            Q(1,Nx) = Vec(Nuse)
         ENDIF
      ENDIF

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

      RETURN
      END