SUBROUTINE DAF_Matrix(N, dx, usys2, X_Values, V, Eig_Values, Eig_Vectors, nd_daf, ndaf, DAFx, dm2, Potential_Name) 
!DAF Matrix
!See ?
USE Numeric_Kinds_Module
IMPLICIT NONE
CHARACTER(LEN=1), PARAMETER:: Jobz='V', UpLo='L'
CHARACTER(LEN=27), PARAMETER:: OutDir='E:\ParkerE\DAF_And_Numerov\'
CHARACTER(LEN=*) Potential_Name
INTEGER, PARAMETER:: lwork=100000, mval=32
INTEGER N, i, j, info, S_work(lwork), IPIV(lwork), nd_daf, ndaf, kmax, NBound
REAL(KIND=WP_Kind) usys2, dx
REAL(KIND=WP_Kind) Eig_Values(N), Eig_Vectors(N,N), X_Values(N), V(N), DAFx(0:nd_daf,0:NDAF), hermit(1000), toe(1000), dm2(n,n)

!Call Get_DAF to find the top row of the DAF matrix
CALL Get_DAF(N,X_Values,DX,Mval,DAFx,Hermit,NDAF,nd_DAF,TOE,usys2,Kmax)

!Since the Matrix is Toeplitz, Symmetric and Banded it is easy to creat the entire matrix
DO i=1,N
    DO j=1,N
        dm2(i,j)=0.d0
        IF(ABS(i-j)<=35)THEN
           dm2(i,j)=DAFx(ABS(i-j),2)
        ENDIF
    ENDDO
ENDDO

!Create the Kinetic Energy Term
dm2=-dm2/usys2

! Add the diagonal potential matrix to the Kinetic Energy to form the Hamiltonian
DO i=1,N
    dm2(i,i)=dm2(i,i)+V(i)
ENDDO

!Diagonalize the Hamiltonian to get the Eigenvalues and Eigenvectors of the system
CALL dsyev(Jobz, Uplo, N, dm2, N, Eig_Values, s_work, lwork, info)

!Write the eigenvalues
OPEN(Unit=131,File=TRIM(OutDIr)//TRIM(Potential_Name)//'_'//"DAF_Eigenvalues.csv")
WRITE(131,*)TRIM(Potential_Name)
WRITE(131,*)"Info=",Info
WRITE(131,*)"DAF Eigenvalues"
WRITE(131,*)"State,  Eig(Hartree),  Eig(cm^-1), Type"
DO I=1,N
   IF(V(N)>=Eig_Values(i))THEN 
      WRITE(131,'(I5,",",ES22.15,",",ES22.15,",",A)')i,Eig_Values(i),Eig_Values(i)*1822.0,"  Bound"
      NBound=I
   ELSE
      WRITE(131,'(I5,",",ES22.15,",",ES22.15,",",A)')i,Eig_Values(i),Eig_Values(i)*1822.0,"  Continuum"
   ENDIF
ENDDO
CLOSE(Unit=131)

!Write the eigenvectors
OPEN(Unit=131,File=TRIM(OutDIr)//TRIM(Potential_Name)//'_'//"DAF_Eigenfunctions.csv")
WRITE(131,*)TRIM(Potential_Name)
WRITE(131,*)"Info=",Info
WRITE(131,*)"DAF Eigenvectors"
DO i=1,N
   WRITE(131,'(201(ES14.5,","))')X_Values(i),(dm2(i,j),j=1,NBound)
ENDDO
CLOSE(Unit=131)

RETURN
ENDSUBROUTINE DAF_Matrix