PROGRAM main
C Simple driver program to check the SLLDEF code
      DOUBLE PRECISION ELAM,EPS,A,B,TOL
      INTEGER K,IFAIL

      WRITE (6,FMT=*) 'Specify the eigenvalue index:'
      READ (5,FMT=*) K
      IF (K.EQ.0) THEN
          WRITE (6,FMT=*) 'Specify positive or negative value of'
          WRITE (6,FMT=*) 'lambda to indicate whether you want to'
          WRITE (6,FMT=*) 'find lambda_{+0} or lambda_{-0}:'
          READ (5,FMT=*) ELAM
      END IF
      EPS = MAX(1.D0,0.1D0*ABS(ELAM))
      WRITE (6,FMT=*) 'Specify tolerance:'
      READ (5,FMT=*) TOL
      A = 0.D0
      B = 1.D0
      WRITE (6,FMT=*) 'Specify endpoints A and B:'
      READ (5,FMT=*) A,B
      IFAIL = 0
      CALL SLLDEF(ELAM,EPS,A,B,K,TOL,IFAIL)
      WRITE (6,FMT=*) 'Eigenvalue approximation found:',ELAM
      STOP
      END

      DOUBLE PRECISION FUNCTION P(X)
      DOUBLE PRECISION X
      P = 1.D0
      RETURN
      END

      DOUBLE PRECISION FUNCTION Q(X)
      DOUBLE PRECISION X
      Q = 0.113785D0
      RETURN
      END

      DOUBLE PRECISION FUNCTION W(X)
      DOUBLE PRECISION X
C      IF (X.LT.0.5D0) THEN
C          W = -1.D0
C      ELSE IF (X.GT.0.5D0) THEN
C          W = 1.D0
C      ELSE
C          W = 0.D0
C      END IF

      W = COS(X)

      RETURN
      END

      SUBROUTINE BCS(C,D,N)
      INTEGER N
      DOUBLE PRECISION C(N,N),D(N,N)

C Neumann boundary conditions:

          C(1,1) = 1.D0
          C(2,1) = 0.D0
          C(1,2) = -1.D0
          C(2,2) = 0.D0
          D(1,1) = 0.D0
          D(2,1) = 1.D0
          D(1,2) = 0.D0
          D(2,2) = 1.D0

      RETURN
      END