      PROGRAM MAIN
      IMPLICIT NONE
C     .. Parameters ..
      INTEGER NZMAX,NFMAX,IMESH,IMLOC,NEV,SIZE,NMESH1
      PARAMETER (NZMAX=5000,NFMAX=NZMAX,IMESH=100000,IMLOC=10,
     & NEV=40,SIZE=2000)
C     ..
C     .. Local Scalars ..
      DOUBLE COMPLEX ANS,ANSBC,ANSLC,ANSRC,ANSTC,
     &               ANSWR1,ANSWR2,DMISS1,DMISS2,EBL,EBR,ELAM,
     &               ELAM1,ELAM2,ETL,ETR,F,ZMIN,FMIN,ROT
      DOUBLE PRECISION A,B,C,CONST,CONST1,CONST2,RELAM,TOL,CN,
     &                 ER,EI,ALOC,BLOC
      INTEGER I,IFAIL,ILOC,IMATCH,IREVCM,NLOC,NMESH,NZACT1,NZACT2,
     &        NZBC,NZLC,NZRC,NZTC,NEIGS,NZ,LSOL,NLOC1,IMAT1
      CHARACTER*1 CYN
C     ..
C     .. Local Arrays ..
      DOUBLE COMPLEX FARRY1(0:NFMAX,1:2),FARRY2(0:NFMAX,1:2),
     &       FB(0:NFMAX,1:2,NEV),FBC1(0:NFMAX,1:2),FL(0:NFMAX,1:2,NEV),
     &       FLC1(0:NFMAX,1:2),FR(0:NFMAX,1:2,NEV),FRC1(0:NFMAX,1:2),
     &       FT(0:NFMAX,1:2,NEV),FTC1(0:NFMAX,1:2),ZARRY1(0:NZMAX),
     &       ZARRY2(0:NZMAX),ZB(0:NZMAX,NEV),ZBC1(0:NZMAX),
     &       ZL(0:NZMAX,NEV),ZLC1(0:NZMAX),ZR(0:NZMAX,NEV),
     &       ZRC1(0:NZMAX),ZT(0:NZMAX,NEV),ZTC1(0:NZMAX),
     &       FBC2(0:NFMAX,1:2),FRC2(0:NFMAX,1:2),FTC2(0:NFMAX,1:2),
     &       FLC2(0:NFMAX,1:2),ZBC2(0:NZMAX),ZLC2(0:NZMAX),
     &       ZRC2(0:NZMAX),ZTC2(0:NZMAX)
      DOUBLE COMPLEX ANSL(NEV),ANSR(NEV),ANST(NEV),ANSB(NEV),
     &       EIGVAL(NEV),ZMAX(2)
      INTEGER NZL(NEV),NZR(NEV),NZT(NEV),NZB(NEV),IEIGS(NEV)
      DOUBLE PRECISION XLOC(0:IMLOC),XMESH(0:IMESH),ARGQHM(2)
      DOUBLE PRECISION YSOL(0:2,1:SIZE),ZSOL(0:2,1:SIZE)
C     ..
C     .. Scalars in Common ..
      DOUBLE PRECISION ALPHA,BETA,NU,EPS
      INTEGER IPROB,PE
      DOUBLE COMPLEX C12,C13,C14,C15,C161,C162,C163,C164
C     ..
C     .. External Subroutines ..
      EXTERNAL BCHILD,BXSIDE,LCHILD,LRCHLD,MKBOX,RCHILD,RFNBOX,SL2BCS,
     &         SL2COF,SLMISS,SLSMBC,SLSMPC,TCHILD,TBCHLD,ALLEVS,GETMSH,
     &         NZERO,STLI02
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC ABS,ATAN,CMPLX,DBLE,AIMAG,NINT
C     ..
C     .. Common blocks ..
      COMMON ALPHA,BETA,NU,IPROB
      COMMON /PROB12/ C12
      COMMON /PROB13/ C13
      COMMON /PROB14/ C14,EPS
      COMMON /PROB15/ C15
      COMMON /PROB16/ C161,C162,C163,C164,PE
C     ..
      NU = 0.D0
      C161 = CMPLX(0.D0,0.D0)
      C162 = CMPLX(0.D0,0.D0)
      C163 = CMPLX(0.D0,0.D0)
      C164 = CMPLX(0.D0,0.D0)
      WRITE (6,FMT=*) 'Specify problem:'
      READ (5,FMT=*) IPROB
      IF (IPROB.LE.0.OR.IPROB.GT.22) STOP
      IF (IPROB.EQ.1) THEN
       BETA = 5.D0
       B = 2.D0*ATAN(1.D0)
       A = -B
      ELSE IF (IPROB.EQ.2) THEN
       A = -20.D0 
       B = 20.D0
      ELSE IF (IPROB.EQ.4) THEN
       A = 0.D0
       B = 10.D0
      ELSE IF (IPROB.EQ.5) THEN
       WRITE (6,FMT=*) 'Specify rotation angle theta:'
       READ (5,FMT=*) BETA
       WRITE (6,FMT=*) 'Specify amplitude parameter C:'
       READ (5,FMT=*) NU
       A = 0.D0
       B = 40.D0
       WRITE (6,FMT=*) 'Specify right endpoint:'
       READ (5,FMT=*) B
      ELSE IF (IPROB.EQ.6) THEN
       A = -1.D0
       B = 1.D0
      ELSE IF (IPROB.EQ.9) THEN
       A = 0.D0
       B = 1.D0
       WRITE (6,FMT=*) 'Specify  amplitude parameter T:'
       READ (5,FMT=*) NU
      ELSE IF (IPROB.EQ.11) THEN
       A = 0.D0
       ALPHA = 1.D0
       BETA = 1.D0
       B = 2.D0*BETA
       NU = 0.D0
       WRITE (6,FMT=*) 'Specify coupling constant:'
       READ (5,FMT=*) NU
      ELSE IF (IPROB.EQ.12) THEN
       WRITE (6,FMT=*) 'Specify Kappa (complex):'
       READ (5,FMT=*) C12
       WRITE (6,FMT=*) 'Specify rotation angle theta:'
       READ (5,FMT=*) NU
       A = 1.D-6
       B = 50.D0
       WRITE (6,FMT=*) 'Specify endpoints:'
       READ (5,FMT=*) A,B
      ELSE IF (IPROB.EQ.13) THEN
       A = 1.D0
       B = 100.D0
       WRITE (6,FMT=*) 'Specify endpoints:'
       READ (5,FMT=*) A,B
       WRITE (6,FMT=*) 'Specify scaling angle:'
       READ (5,FMT=*) NU
       WRITE (6,FMT=*) 'Specify power BETA:'
       READ (5,FMT=*) BETA
      ELSE IF (IPROB.EQ.14) THEN
       WRITE (6,FMT=*) 'Specify Kappa (complex):'
       READ (5,FMT=*) C14
       WRITE (6,FMT=*) 'Specify P:'
       READ (5,FMT=*) PE
       WRITE (6,FMT=*) 'Specify rotation angle theta:'
       READ (5,FMT=*) NU
       A = 1.D-3
       B = 50.D0
       WRITE (6,FMT=*) 'Specify endpoints:'
       READ (5,FMT=*) A,B
       EPS = A
C Change to the left-endpoint for the T-variable, defined
C by X = EPS + T*EXP(I*NU):
       A = 0.D0
      ELSE IF (IPROB.EQ.15) THEN
       WRITE (6,FMT=*) 'Specify Kappa (complex):'
       READ (5,FMT=*) C15
       WRITE (6,FMT=*) 'Specify rotation angle theta:'
       READ (5,FMT=*) NU
       A = -10.D0
       B = 10.D0
       WRITE (6,FMT=*) 'Specify endpoints:'
       READ (5,FMT=*) A,B       
      ELSE IF (IPROB.EQ.16) THEN
       WRITE (6,FMT=*) 'Specify endpoints:'
       READ (5,FMT=*) A,B       
       WRITE (6,FMT=*) 'Specify rotation angle theta:'
       READ (5,FMT=*) NU
       WRITE (6,FMT=*) 'Specify Kappa (complex):'
       READ (5,FMT=*) C162
       WRITE (6,FMT=*) 'Specify coupling constant (complex):'
       READ (5,FMT=*) C161
       WRITE (6,FMT=*) 'Specify power P:'
       READ (5,FMT=*) PE
      ELSE IF (IPROB.EQ.17) THEN
       WRITE (6,FMT=*) 'Specify endpoints:'
       READ (5,FMT=*) A,B    
       WRITE (6,FMT=*) 'Specify C, Kappa (both complex):'
       READ (5,FMT=*) C161,C162
       WRITE (6,FMT=*) 'Specify rotation angle:'
       READ (5,FMT=*) NU
       C161 = CMPLX(COS(NU),SIN(NU))*C161   
       C162 = CMPLX(COS(NU),SIN(NU))*C162   
       PE = 1
      ELSE IF (IPROB.EQ.19) THEN
       WRITE (6,FMT=*) 'Specify endpoints:'
       READ (5,FMT=*) A,B    
       WRITE (6,FMT=*) 'Specify C, Kappa (both complex):'
       READ (5,FMT=*) C161,C162
       WRITE (6,FMT=*) 'Specify rotation angle:'
       READ (5,FMT=*) NU
       ROT = CMPLX(COS(NU),SIN(NU))
       C161 = C161/(ROT**2)
       C162 = C162*ROT
      ELSE IF (IPROB.EQ.20) THEN
       WRITE (6,FMT=*) 'Specify endpoints:'
       READ (5,FMT=*) A,B    
       WRITE (6,FMT=*) 'Specify C, Kappa, Alfa, Beta:'
       READ (5,FMT=*) C161,C162,C163,C164
       WRITE (6,FMT=*) 'Specify rotation angle:'
       READ (5,FMT=*) NU
       ROT = CMPLX(COS(NU),SIN(NU))
       C161 = C161*ROT
       C162 = C162*ROT
       C163 = C163*ROT
       C164 = C164*(ROT**2)
      ELSE IF (IPROB.EQ.21) THEN
       WRITE (6,FMT=*) 'Specify endpoints:'
       READ (5,FMT=*) A,B    
       WRITE (6,FMT=*) 'Specify C, Kappa, Alfa, Beta:'
       READ (5,FMT=*) C161,C162,C163,C164
      ELSE IF (IPROB.EQ.22) THEN
       WRITE (6,FMT=*) 'Specify endpoints:'
       READ (5,FMT=*) A,B    
       WRITE (6,FMT=*) 'Specify C, Kappa, Alfa, Beta:'
       READ (5,FMT=*) C161,C162,C163,C164
      ELSE
       A = 0.D0
       B = 1.D0
       WRITE (6,FMT=*) 'Specify endpoints:'
       READ (5,FMT=*) A,B
      END IF

      WRITE (6,FMT=*) 'Specify meshing tolerance:'
      READ (5,FMT=*) TOL
      NMESH = IMESH
      IFAIL = 0
      ELAM = CMPLX(-5.D0,1.D0)
      ALOC = A
      BLOC = B


C      IF (IPROB.EQ.16) THEN
C          ALOC = A/(1.D0+ABS(A))
C          BLOC = B/(1.D0+ABS(B))
C          NMESH = MIN(IMESH,NINT((B-BLOC)*TOL**(-0.25D0)))     
C          DO 2003 I = 0,NMESH
C             XMESH(I) = ALOC + (BLOC-ALOC)*DBLE(I)/DBLE(NMESH)
C             XMESH(I) = XMESH(I)/(1.D0-ABS(XMESH(I)))
C2003      CONTINUE
C          IMATCH = 0
C          GOTO 2004
C      END IF

      IF (IPROB.EQ.12.OR.IPROB.EQ.14) THEN
         ALOC = A
         BLOC = 1.D0+MAX(0.D0,A)
         WRITE (6,FMT=*) '1st part of meshing over',ALOC,BLOC
      END IF 
   
      CALL GETMSH(ALOC,BLOC,XMESH,NMESH,IMATCH,ELAM,SL2COF,SL2BCS,TOL,
     &              IFAIL)
      IF (IFAIL.NE.0) THEN
         WRITE (6,FMT=*) 'Meshing failed, IFAIL =',IFAIL
         STOP
      END IF
      WRITE (6,FMT=*) 'Mesh has',NMESH,'intervals'
      WRITE (7,FMT=*) 'Here is the mesh:'
      DO 173 I = 0,NMESH
         WRITE (7,FMT=*) I,XMESH(I)
173   CONTINUE

      IF (IPROB.EQ.12.OR.IPROB.EQ.14) THEN
         IMATCH = NMESH  
      END IF      

      IF (IPROB.EQ.12.OR.IPROB.EQ.14) THEN
      ALOC = XMESH(IMATCH)
      BLOC = XMESH(IMATCH)
      NLOC1 = IMESH
      CALL GETMSH(A,BLOC,XMESH,NLOC1,IMAT1,ELAM,SL2COF,SL2BCS,TOL,IFAIL)
      IF (IFAIL.NE.0) THEN
         WRITE (6,FMT=*) 'Meshing failed II, IFAIL =',IFAIL
         STOP
      END IF

      IMATCH = NLOC1

      NMESH = NINT((B-BLOC)*TOL**(-0.25D0))
      IF (NMESH+IMATCH.GT.IMESH) NMESH = IMESH-IMATCH
      WRITE (6,FMT=*) 'IMATCH,NMESH,IMESH:',IMATCH,NMESH,IMESH
      DO 33 I = 1,NMESH
       XMESH(IMATCH+I) = BLOC+I*(B-BLOC)/NMESH
33    CONTINUE
      NMESH = IMATCH+NMESH
      END IF

      IF (IPROB.EQ.6) THEN
      NMESH = 10
      IMATCH = 5
      DO 32 I = 0,NMESH
         XMESH(I) = 2.D0*DBLE(I)/DBLE(NMESH)-1.D0
32    CONTINUE
      END IF

      IF (IPROB.EQ.11) THEN
      NMESH = 10
      IMATCH = 5
      DO 321 I = 0,NMESH
         XMESH(I) = 2.D0*DBLE(I)/DBLE(NMESH)
321    CONTINUE
      END IF

      IF (IPROB.EQ.21.OR.IPROB.EQ.20.OR.IPROB.EQ.19.OR.
     &    IPROB.EQ.17.OR.IPROB.EQ.5) THEN
      NMESH1 = NINT(B)
      WRITE (6,FMT=*) 'Specify NMESH:'
      READ (5,FMT=*) NMESH1
      IF (NMESH1.LE.0) GOTO 2004
      NMESH = NMESH1
      IMATCH = 1
C      IF (NU.NE.0.D0) THEN
C          IMATCH = 1
C      ELSE
C          IMATCH = NMESH-1
C      END IF
      DO 322 I = 0,NMESH
         XMESH(I)=A+I*(B-A)/NMESH
322   CONTINUE
      END IF

2004  CONTINUE

      NLOC = 10
      DO 30 I = 0,NLOC
          XLOC(I) = (DBLE(NLOC-I)*A+DBLE(I)*B)/DBLE(NLOC)
   30 CONTINUE
      ILOC = NINT(DBLE(NLOC)/2.D0)

      IF (IPROB.EQ.4) THEN
        ILOC = 0
      END IF

      IF (IPROB.NE.4.AND.IPROB.NE.6) THEN
        WRITE (6,FMT=*) 'Do you want to do a half-plane count (Y/N)?'
        READ (5,FMT=9000) CYN
      ELSE
        CYN = 'N'
      END IF
9000  FORMAT(A)

      IF (CYN.EQ.'y'.OR.CYN.EQ.'Y') THEN

      WRITE (6,FMT=*) 'Counting all eigenvalues with real part less',
     &                ' than RELAM; specify RELAM:'
      READ (5,FMT=*) RELAM

      IFAIL = 0
C
C NZERO finds all eigenvalues with real part less than RELAM. In order to
C use this routine you must be able to find a second eigenvalue problem 
C for which you know (a) that there are no eigenvalues with real part
C less than or equal to RELAM; (b) the **total order** of the boundary 
C conditions is the same as for the original problem in which you are
C interestes. For example, if your problem has boundary conditions
C
C   y'(0) - y(0) = 0,   y(1) = 0,
C
C then the total order of the boundary conditions is 1 + 0 = 0 (the
C sums of the highest orders of the derivatives appearing). In this case
C an appropriate second eigenproblem might use the simpler boundary
C conditions 
C
C   y'(0) = 0, y(1) = 0,
C
C for which it is easier to find a problem all of whose eigenvalues are
C to the right of RELAM.
C
C This second problem must be specified by the routines SLSMPC and SLSMBC
C (coefficients and boundary conditions respectively) just as the original
C problem in which you are really interested is specified by SL2COF and
C SL2BCS.
C

      CALL NZERO(XMESH,NMESH,IMATCH,RELAM,SL2COF,SL2BCS,SLSMPC,
     &           SLSMBC,NZ,ZMIN,FMIN,ZMAX,ARGQHM,IFAIL)

      IF (IFAIL.NE.0) THEN
          WRITE (6,FMT=*) 'Failure in call to NZERO, IFAIL=',IFAIL
          STOP
      END IF

      WRITE (6,FMT=*) 'No. of eigenvalues with real part less than',
     & RELAM,' is',NZ

      END IF

      WRITE (6,FMT=*) 'Do you want to find all the eigenvalues in',
     & ' a box (Y/N)?'
      READ (5,FMT=9000) CYN
      IF (CYN.NE.'Y'.AND.CYN.NE.'y') STOP

C Now decide how many eigenvalues are in a given box.
      WRITE (6,FMT=*) 'Specify bottom right corner of box:'
      READ (5,FMT=*) EBR
      WRITE (6,FMT=*) 'Specify top left corner of box:'
      READ (5,FMT=*) ETL
      ETR = CMPLX(0.D0,0.D0)
      EBL = CMPLX(0.D0,0.D0)
      ETR = CMPLX(DBLE(EBR),AIMAG(ETL))
      EBL = CMPLX(DBLE(ETL),AIMAG(EBR))

      WRITE (6,FMT=*) 'Top right corner:',ETR
      WRITE (6,FMT=*) 'Bottom left corner:',EBL

      IFAIL = 0
C
C The routine MKBOX integrates round a box in the complex plane with corners
C EBR (bottom right), ETR (top right), ETL (top left) and EBL (bottom left),
C thereby counting the eigenvalues in the box. It must be called before any
C attempt to find the eigenvalues in the box.

      CALL MKBOX(EBR,ETR,ETL,EBL,ZL,NZMAX,NZL,ZR,NZMAX,NZR,ZT,NZMAX,NZT,
     &           ZB,NZMAX,NZB,FL,NFMAX,FR,NFMAX,FT,NFMAX,FB,NFMAX,ANSL,
     &           ANSR,ANST,ANSB,ANS,XMESH,NMESH,IMATCH,SL2COF,SL2BCS,
     &           IFAIL)
      NEIGS = NINT(AIMAG(ANS))
      WRITE (6,FMT=*) 'IFAIL was',IFAIL
      WRITE (6,FMT=*) 'Number of eigenvalues was',NEIGS
      WRITE (6,FMT=*) 'ANS was',ANS
      WRITE (6,FMT=*) 'Numbers NZR,NZT,NZL,NZB:',NZR(1),NZT(1),NZL(1),
     &                 NZB(1)

      IF (NEIGS.LE.NEV) THEN
      WRITE (6,FMT=*) 'Computing all eigenvalues:'
      WRITE (6,FMT=*) 'Specify tolerance:'
      TOL = 1.D-7
      READ (5,FMT=*) TOL
      IFAIL = 0
      DO 77 I = 1,NEIGS
         EIGVAL(I) = CMPLX(0.D0,0.D0)
77    CONTINUE
      IEIGS(1) = NEIGS
C
C The routine ALLEVS can be called after a call to MKBOX to find all the
C eigenvalues in the box specified. The storage arrays supplied must be
C big enough for the number of eigenvalues in the box, hence the 
C IF (NEIGS.LE.NEV) statement above. TOL controls the rootfinding
C accuracy and EIGVAL(1)...EIGVAL(NEIGS) will contain the eigenvalues
C on exit.
C

      CALL ALLEVS(ZL,NZMAX,NZL,ZR,NZMAX,NZR,ZT,NZMAX,NZT,ZB,NZMAX,NZB,
     &            FL,NFMAX,FR,NFMAX,FT,NFMAX,FB,NFMAX,ANSL,ANSR,ANST,
     &            ANSB,ZLC1,NZMAX,ZRC1,NZMAX,ZTC1,NZMAX,ZBC1,NZMAX,
     &            FLC1,NFMAX,FRC1,NFMAX,FTC1,NFMAX,FBC1,NFMAX,
     &            ZLC2,NZMAX,ZRC2,NZMAX,ZTC2,NZMAX,ZBC2,NZMAX,
     &            FLC2,NFMAX,FRC2,NFMAX,FTC2,NFMAX,FBC2,NFMAX,
     &            XMESH,NMESH,IMATCH,SL2COF,SL2BCS,TOL,
     &            EIGVAL,NEIGS,IEIGS,IFAIL)

      WRITE (6,FMT=*) 'Completed problem',IPROB

      IF (IPROB.NE.5) THEN
      WRITE (6,FMT=*) 'Here are the eigenvalues:'
      DO 17 I = 1,NEIGS
      ER = DBLE(EIGVAL(I))
      EI = AIMAG(EIGVAL(I))
      WRITE (6,FMT=*) EIGVAL(I)
C ,SQRT(ER**2+EI**2),
C     &                ATAN2(EI,ER)
      WRITE (16,FMT=*) EIGVAL(I)

17    CONTINUE      
      WRITE (6,FMT=*) 'IFAIL was',IFAIL
      ELSE
      WRITE (6,FMT=*) 'Here are the eigenvalues (values of z^2):'
      DO 19 I = 1,NEIGS
      WRITE (6,FMT=*) EIGVAL(I)/CMPLX(COS(2.D0*BETA),
     &                                 SIN(2.D0*BETA))
      WRITE (16,FMT=*) EIGVAL(I)/CMPLX(COS(2.D0*BETA),
     &                                 SIN(2.D0*BETA))
19    CONTINUE   
      WRITE (6,FMT=*) 'Here are the resonances (values of z):'
      DO 191 I = 1,NEIGS
      WRITE (6,FMT=*) (EIGVAL(I)/CMPLX(COS(2.D0*BETA),
     &                                 SIN(2.D0*BETA)))**0.5D0
      WRITE (16,FMT=*) (EIGVAL(I)/CMPLX(COS(2.D0*BETA),
     &                                 SIN(2.D0*BETA)))**0.5D0
191   CONTINUE    
      WRITE (16,FMT=*) 'Eigenvalues of rotated problem:'
      DO 192 I = 1,NEIGS
      WRITE (16,FMT=*) EIGVAL(I)
192   CONTINUE    


      WRITE (6,FMT=*) 'IFAIL was',IFAIL
      END IF
      IF (IPROB.EQ.12) WRITE (6,FMT=*) 'Matchpoint:',XMESH(IMATCH)
C      WRITE (6,FMT=*) 'Specify which eigenfunction to compute:'
C      READ (5,FMT=*) I
C      IF (I.LT.1.OR.I.GT.NEIGS) STOP 'Eigenvalue choice out of range'
C      WRITE (6,FMT=*) 'Tolerance for eigenfunction calculation:',TOL
C      WRITE (6,FMT=*) 'Specify matching point:'
C      READ (5,FMT=*) C
C      IF (C.LT.A.OR.C.GT.B) C = XMESH(IMATCH)
C      
C      CALL STLIO2(EIGVAL(I),A,B,C,SIZE,TOL,SL2BCS,
C     &            SL2COF,YSOL,ZSOL,LSOL,CN)
C
C      DO 10 I = 1,LSOL
C        WRITE (27,FMT=*) YSOL(0,I),YSOL(1,I),YSOL(2,I)
C        WRITE (28,FMT=*) ZSOL(0,I),ZSOL(1,I),ZSOL(2,I)
C10    CONTINUE
C      WRITE (6,FMT=*) 'Condition number was',CN

      END IF

      STOP
      END

C ------------------------------------------------------------------------------



      SUBROUTINE SL2COF(X,ELAM,P,Q1,Q0)
      IMPLICIT NONE
C 
C This subroutine specifies the coefficients in the differential equation. The
C differential equation is
C
C   -(P(x)y')' + Q1(x)y' + Q0(x,lambda)y = 0.
C
C X denotes x, ELAM denotes lambda. The user must write this routine to specify
C the values of the coefficients which correspond to a given X and ELAM.
C
C
C     .. Scalar Arguments ..
      DOUBLE COMPLEX ELAM,Q0,Q1
      DOUBLE PRECISION P,X
C     ..
C     .. Local Scalars ..
      DOUBLE PRECISION Q,S,PI
      DOUBLE COMPLEX C,Z,QL0,QL1,DQL0,D2QL0,R,DR,D2R
C     ..
C     .. External Functions ..
      DOUBLE COMPLEX ZED,QBASE,XI
      DOUBLE PRECISION X01AAF
      EXTERNAL ZED,QBASE,XI
      EXTERNAL GETRT
      EXTERNAL X01AAF
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC CMPLX,COS,SIN
C     ..

C     .. Scalars in Common ..
      DOUBLE PRECISION ALPHA,BETA,NU,EPS
      INTEGER IPROB,PE
      DOUBLE COMPLEX C12,C13,C14,C15,C161,C162,C163,C164
      DOUBLE COMPLEX CONST(17)
C     ..
C     .. Common blocks ..
      COMMON ALPHA,BETA,NU,IPROB
      COMMON /PROB12/ C12
      COMMON /PROB13/ C13
      COMMON /PROB14/ C14,EPS
      COMMON /PROB15/ C15
      COMMON /PROB16/ C161,C162,C163,C164,PE
C     ..
      GOTO (10,20,30,40,50,60,70,80,90,100,110,
     &      120,130,140,150,160,170,180,190,200,
     &      210,220) IPROB

C The first problem is the Coffey-Evans equation with a complex shifted potential.
10    P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Q = -2.D0*BETA*COS(2.D0*X) + (BETA*SIN(2.D0*X))**2
      Q0 = CMPLX(Q,1.D0) - ELAM
      RETURN

C The second problem is Brian Davies rotated complex harmonic oscillator.
20    P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Q0 = CMPLX(X**2,X**3+(X**2)*SIN(X))- ELAM
      RETURN

C A Schroedinger problem on a compact interval
30    P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      PI = X01AAF(0.D0)
C      CONST(1) = (3.0015999800525606E-04,1.1535248719155788E-02)
C      CONST(2) = (-2.7713517192751169E-04,-1.2606701347976923E-03)
C      CONST(3) = (-5.1086430175928399E-05,-2.3820089700166136E-04)
C      CONST(4) = (3.2295236451318488E-05,1.5909824287518859E-04)
C      CONST(5) = (-2.2328862542053685E-05,-1.1323267972329631E-04)
C      CONST(6) = (1.6455027434858494E-05,8.4203114965930581E-05)
C      CONST(7) = (-1.2584612704813480E-05,-6.4921827288344502E-05)
C      CONST(8) = (9.9314192993915640E-06,5.1531929784687236E-05)
C      CONST(9) = (-8.0341878856415860E-06,-4.1876395698636770E-05)
C      CONST(10) = (3.2607156783342361E-02,-8.2908503711223602E-02)
C      CONST(11) = (-7.8643383458256721E-03,-3.3953104168176651E-02)
C      CONST(12) = (7.8567313030362129E-03,5.1450252532958984E-02)
C      CONST(13) = (-8.2074701786041260E-03,-4.6013932675123215E-02)
C      CONST(14) = (8.1627769395709038E-03,4.4367637485265732E-02)
C      CONST(15) = (-8.1037599593400955E-03,-4.3386802077293396E-02)
C      CONST(16) = (8.0157453194260597E-03,4.2896773666143417E-02)
C      CONST(17) = (-7.9882368445396423E-03,-4.2610663920640945E-02)
      Q0 = CMPLX(1.D0,-1.D0+10.D0*COS(PI*X)) - ELAM
C      Q0 = CONST(1) + CONST(2)*SIN(X)+CONST(3)*COS(X)
C     &   + CONST(4)*SIN(2.D0*X)+CONST(5)*COS(2.D0*X)
C     &   + CONST(6)*SIN(3.D0*X)+CONST(7)*COS(3.D0*X)
C     &   + CONST(8)*SIN(4.D0*X)+CONST(8)*COS(4.D0*X)
C     &   + CONST(9)*SIN(5.D0*X)+CONST(10)*COS(5.D0*X)
C     &   + CONST(11)*SIN(6.D0*X)+CONST(12)*COS(6.D0*X)
C     &   + CONST(13)*SIN(7.D0*X)+CONST(14)*COS(7.D0*X)
C     &   + CONST(15)*SIN(8.D0*X)+CONST(16)*COS(8.D0*X)
C     &   + CONST(17)*SIN(9.D0*X) - ELAM

      RETURN

C This is -y'' + exp(-x^2)y = \lambda y
40    P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Q0 = CMPLX(-EXP(-X*X),0.D0) - ELAM
      RETURN      

C This is a resonance problem for Malcolm Brown
50    P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      C = CMPLX(COS(BETA),SIN(BETA))
      Z = C*CMPLX(X,0.D0)
C      Q0 = CMPLX(-16.D0*NU,0.D0)*EXP(-Z/4.D0)*COS(Z)
C      Q0 = Z*Z*EXP(-Z*Z/5.D0)
      Q0 = 5.D0*Z*EXP(-Z)
C      Q0 = Z*Z*Z*EXP(-Z*Z*Z/100.D0)
C      Q0 = -16.D0*Z*Z*EXP(-Z)
C      IF (X.GE.40.D0) Q0 = CMPLX(0.D0,0.D0)
C      Q0 = Z*Z*EXP(-NU*Z*Z)
C      Q0 = 1.D0/(CMPLX(1.D0,0.D0)+Z)**2
C      Q0 = NU*(Z-CMPLX(1.D0,0.D0))/((Z+CMPLX(1.D0,0.D0))**4)
C      Q0 = CMPLX(0.25D0,0.D0)*EXP(-2.0D0*Z)
      Q0 = C*C*Q0 - ELAM
      RETURN
C This is a left-definite problem
60    P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Q0 = -ELAM*SIGN(1.D0,X)
C      WRITE (6,FMT=*) 'X,P,Q0:',X,P,Q0
      RETURN
C
70    P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      C = CMPLX(1.D0,1.D0)
      Q0 = C*C-ELAM*EXP(-3.D0*X*ABS(DBLE(C)))
C      WRITE (6,FMT=*) 'X,P,Q0:',X,P,Q0
      RETURN

80    P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Q0 = CMPLX(0.D0,X**(-0.75D0)) - ELAM
 
      RETURN 

90    P = 1.D0
      PI = 4.D0*ATAN(1.D0)
      Q1 = CMPLX(0.D0,-PI)
      Q0 = (PI*PI/4.D0)*CMPLX(COS(2.D0*PI*X),-SIN(2.D0*PI*X))
     &     *QBASE(ZED(X),ELAM,NU)
      RETURN

100   P = 1.D0
      Q1 = CMPLX(-X,2.D0)/CMPLX(1.D0,X)
      Q0 = -ELAM
      RETURN

110   P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      IF (X.LT.BETA) THEN
         Q = ALPHA + NU
      ELSE
         Q = 0.D0
      END IF
      Q0 = CMPLX(Q,0.D0)-ELAM
C      WRITE (11,FMT=*) X,Q0,ALPHA,BETA,NU
C      WRITE (11,FMT=*)
      RETURN

120   P = 1.D0
      C = CMPLX(COS(NU),SIN(NU))
      Q1 = CMPLX(0.D0,0.D0)
      Z = XI(C12,C*X)
C      Q0 = C*C*( 1.D0/Z - C12**2 + (ELAM*ELAM+ELAM)/(Z*Z) )
      QL0 = C12**2 - 1.D0/(1.D0+Z)+0.25D0/((C12*(1.D0+Z))**2)
      DQL0 = 1.D0/(1.D0+Z)**2 - 0.5D0/((C12**2)*((1.D0+Z)**3))
      D2QL0 = -2.D0/((1.D0+Z)**3) - 1.5D0/((C12**2)*((1.D0+Z)**4))
      QL1 = 1.D0/(Z*(Z+1.D0))-ELAM*(ELAM+1.D0)/(Z*Z) 
     &    - 0.25D0/((C12*(Z+1.D0))**2)
      Q0 = 1.D0+(5.D0/16.D0)*((DQL0/QL0)**2)/QL0 
     &   - 0.25D0*D2QL0/(QL0*QL0) + QL1/QL0
      Q0 = -Q0
     
      RETURN

130   P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Q0 = -CMPLX(COS(NU),SIN(NU))/(X**BETA)-ELAM
      RETURN

140   P = 1.D0
      C = CMPLX(COS(NU),SIN(NU))
      Q1 = CMPLX(0.D0,0.D0)
      Z = EPS+C*X
      Q0 = C*C*( 1.D0/Z**PE -C14**2 + (ELAM*ELAM+ELAM)/(Z*Z) )
C      WRITE (7,FMT=*) 'EPS,Kappa,NU:',EPS,C14,NU
      RETURN

150   P = 1.D0
      C = CMPLX(COS(NU),SIN(NU))
      Q1 = CMPLX(0.D0,0.D0)
      QL0 = EXP(X)*C
      Q0 = QL0**2 
      Q0 = QL0 - (C15**2)*Q0 + CMPLX(0.25D0,0.D0) 
     &   + ELAM*ELAM+ELAM
      RETURN

160   P = 1.D0
      C = CMPLX(COS(NU),SIN(NU))
      Q1 = CMPLX(0.D0,0.D0)
      QL0 = EXP(X)*C
      Q0 = QL0**2 
      Q0 = CMPLX(COS((2.D0-PE)*NU),SIN((2.D0-PE)*NU))
     &                          *EXP((2.D0-PE)*X)*C161
     &   - (C162**2)*Q0 + CMPLX(0.25D0,0.D0) 
     &   + ELAM*ELAM+ELAM
      RETURN

170   P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Z = CMPLX(X,0.D0)
      Q0 = C161/Z - C162**2 + (ELAM*ELAM+ELAM)/(Z*Z) 

      RETURN

180   P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Q0 = CMPLX(1.D0,1.D0)*X**3-ELAM

      RETURN

190   P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Q0 = C161/(X**4) -C162**2 + (ELAM+ELAM**2)/(X**2)
      RETURN

200   P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Z = CMPLX(X,0.D0)
      Q0 = C161/(Z*(1.D0+C163*Z)*(1.D0+C164*Z*Z)) - C162**2 
     &      + (ELAM*ELAM+ELAM)/(Z*Z)
C      IF (ABS(C163).NE.0.D0.OR.ABS(C164).NE.0.D0) THEN
C
C      WRITE (17,FMT=*) 'C161, C162, C163, C164, Z, ELAM:',
C     & C161,C162,C163,C164,Z,ELAM
C
C      END IF


      RETURN

210   P = 1.D0
      CALL GETRT(X,C161,C162,C163,C164,ELAM,R,DR,D2R)
      Q1 = D2R/DR
      Q0 = (DR**2)*(C161/(R*(1.D0+C163*R)*(1.D0+C164*R*R)) - C162**2 
     &      + (ELAM*ELAM+ELAM)/(R*R))

      RETURN

220   P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Z = CMPLX(X,0.D0)
      Q0 = C161/(Z*(1.D0+C163*Z)*(1.D0+C164*Z*Z)) - C162**2 
     &      + (ELAM*ELAM+ELAM)/(Z*Z)

      END

      SUBROUTINE GETRT(X,C161,C162,C163,C164,ELAM,R,DR,D2R)
      IMPLICIT NONE
      DOUBLE PRECISION X,SCAL,CT,DENOM
      DOUBLE COMPLEX C161,C162,C163,C164,ELAM,R,DR,D2R

      SCAL = ABS(C161*C164*0.25D0)/(ABS(ELAM)**2+1.D0)
      CT = ABS(C164)*X
      DENOM = 1.D0+0.5D0*(CT**3)
      R = CMPLX(SCAL*X/DENOM,X)
      DR = CMPLX(SCAL*(1.D0-CT**3)/(DENOM**2),1.D0)
      D2R = CMPLX( SCAL*(-3.D0*CT**3/(X*(DENOM**2))
     &            -3.D0*(1.D0-CT**3)*(CT**3)/(X*(DENOM**3))),0.D0)
      RETURN
      END

      DOUBLE COMPLEX FUNCTION QBASE(Z,ELAM,NU)      
      DOUBLE COMPLEX Z,ELAM
      DOUBLE PRECISION NU
      QBASE = NU*(CMPLX(2.5D0,0.D0)-Z)*
     &   EXP(-Z)/((Z-CMPLX(0.5D0,0.D0))*EXP(-Z)-ELAM)
     &   + ELAM
      RETURN 



      END

      DOUBLE COMPLEX FUNCTION ZED(X)
      DOUBLE PRECISION X,PI
      PI = 4.D0*ATAN(1.D0)
      ZED = CMPLX(0.5D0*(1.D0-COS(PI*X)),
     &             0.5D0*SIN(PI*X))
      RETURN
      END

C -----------------------------------------------------------

      DOUBLE COMPLEX FUNCTION XI(K,T)
      IMPLICIT NONE
      DOUBLE COMPLEX T
      DOUBLE COMPLEX K
      DOUBLE COMPLEX X,XN

      X = T/K

10    XN = (T+0.5D0*LOG(1.D0+X)/K)/K
      IF (ABS(X-XN).GT.1.D-15*MAX(1.D0,ABS(T))) THEN
          X = XN
          GOTO 10
      END IF
      XI = XN
C      WRITE (17,FMT=*) 'T,X',T,X
      RETURN
      END

C -----------------------------------------------------------

      SUBROUTINE SL2BCS(IEND,X,ELAM,Y,PDY)
      IMPLICIT NONE
C     .. Scalar Arguments ..
      DOUBLE COMPLEX ELAM,Y,PDY
      DOUBLE PRECISION X
      INTEGER IEND
C     ..
C     .. Scalars in Common ..
      DOUBLE PRECISION ALPHA,BETA,NU,EPS
      INTEGER IPROB,PE
      DOUBLE COMPLEX C12,C13,C14,C15,C161,C162,C163,C164
C     ..
C     .. Local Scalars ..
      DOUBLE COMPLEX Z,ALF,A1,A2,A3,A4,QLOC,SQLOC,DQLOC,QL0,
     &               DQL0,Q,DQ,D2Q,DS0,DS1,DS2,R,DR,D2R
      LOGICAL FRSTME
C     ..
C     .. Common blocks ..
      COMMON ALPHA,BETA,NU,IPROB
      COMMON /PROB12/ C12
      COMMON /PROB14/ C14,EPS
      COMMON /PROB15/ C15
      COMMON /PROB16/ C161,C162,C163,C164,PE
C     ..
C     .. External Functions ..
      DOUBLE COMPLEX MYSQRT,XI
      EXTERNAL MYSQRT,XI
C     .. Local Scalars ..
      DOUBLE PRECISION THETA,P,PI,XBIG
      DOUBLE COMPLEX IOMEG1,IOMEG2,E1,E2,Q1,U,DU,ROT
      DOUBLE COMPLEX D,G,LAMBDM,LAMBDP,LAMBDA,Y22B,IX
      DOUBLE COMPLEX B1,B2,B3,IM,KAPPA,CEE,C
      DATA FRSTME /.TRUE./
      SAVE XBIG,FRSTME



      GOTO (10,20,30,40,50,60,70,80,90,100,110,
     &      120,130,140,150,160,170,180,190,200,
     &      210,220) IPROB

C Dirichlet BCs for problem 1:
10    Y = CMPLX(0.D0,0.D0)
      PDY = CMPLX(1.D0,0.D0)

      RETURN
C Dirichlet BCs for problem 2:
20    Y = CMPLX(0.D0,0.D0)
      PDY = CMPLX(1.D0,0.D0)

      RETURN

C Neumann BCs for problem 3:
30    Y = CMPLX(1.D0,0.D0)
      PDY = CMPLX(0.D0,0.D0)

      RETURN
C
C Neumann BCs at left hand end, weird lambda-dependent BCs at right hand
C end for problem 4:
40    IF (IEND.EQ.0) THEN
          Y = CMPLX(1.D0,0.D0)
          PDY = CMPLX(0.D0,0.D0)
      ELSE
          Y = CMPLX(1.D0,0.D0)
          PDY = CMPLX(0.D0,1.D0)*MYSQRT(ELAM)
C          WRITE (10,FMT=*) 'LAMBDA was',ELAM
C          WRITE (10,FMT=*) 'MYSQRT(LAMBDA) was',MYSQRT(ELAM)
      END IF

      RETURN

C Various boundary conditions for problem 5, the resonance problem:
50    IF (IEND.EQ.0) THEN
          Y = CMPLX(1.D0,0.D0)
          PDY = CMPLX(0.D0,0.D0)
C          Y = CMPLX(0.D0,0.D0)
C          PDY = CMPLX(1.D0,0.D0)

      ELSE
C          CALL SL2COF(X,ELAM,P,Q1,Q)
C          Y = CMPLX(1.D0,0.D0)
C          PDY = -SQRT(Q)
C          IF (DBLE(PDY).GT.0.D0) PDY = -PDY
          Y = CMPLX(0.D0,0.D0)
          PDY = CMPLX(1.D0,0.D0)
      END IF

      RETURN

60    Y = CMPLX(0.D0,0.D0)
      PDY = CMPLX(1.D0,0.D0)
      RETURN

70    IF (IEND.EQ.0) THEN
       Y = CMPLX(0.D0,0.D0)
       PDY = CMPLX(1.D0,0.D0)
      ELSE
       Y = CMPLX(1.D0,0.D0)
       PDY = CMPLX(-1.D0,0.D0)
      END IF
      RETURN

80    IF (IEND.EQ.0) THEN
       Y = CMPLX(0.D0,0.D0)
       PDY = CMPLX(1.D0,0.D0)
      ELSE
       Y = CMPLX(1.D0,0.D0)
       PDY = -SQRT(ELAM) 
       IF (DBLE(PDY).GT.0.D0) PDY = -PDY
      END IF
      RETURN

90    IF (IEND.EQ.0) THEN
         U = CMPLX(-0.5D0,0.D0)
        DU = CMPLX(1.5D0,0.D0)
      ELSE
         U = CMPLX(0.5D0,0.D0)
        DU = CMPLX(0.5D0,0.D0)
      END IF
      PDY = DU
      PI = 4.D0*ATAN(1.D0)
      IF (IEND.EQ.0) THEN
        Y = CMPLX(0.D0,-2.D0/PI)*U
      ELSE
        Y = CMPLX(0.D0,2.D0/PI)*U
      END IF

      RETURN

100   IF (IEND.EQ.0) THEN
          Y = CMPLX(0.D0,0.D0)
          PDY = CMPLX(1.D0,0.D0)
      ELSE
           Y = CMPLX(1.D0,0.D0)
          PDY = CMPLX(1.D0,1.D0)
      END IF
      RETURN

110   G = SQRT(ELAM-ALPHA)/SQRT(ELAM)
      D = 2.D0*COS(SQRT(ELAM-ALPHA)*BETA)*COS(SQRT(ELAM)*BETA)
     & - (G+1.D0/G)*SIN(SQRT(ELAM-ALPHA)*BETA)*SIN(SQRT(ELAM)*BETA)
      LAMBDP = (D + SQRT(D*D-4.D0))/2.D0 
      LAMBDM = (D - SQRT(D*D-4.D0))/2.D0 

      Y = -SIN(SQRT(ELAM-ALPHA)*BETA)*
     & COS(SQRT(ELAM)*BETA)/SQRT(ELAM-ALPHA)
     & -SIN(SQRT(ELAM)*BETA)*COS(SQRT(ELAM-ALPHA)*BETA)/SQRT(ELAM)

      Y22B = COS(SQRT(ELAM)*BETA)*COS(SQRT(ELAM-ALPHA)*BETA)-G*
     & SIN(SQRT(ELAM-ALPHA)*BETA)*SIN(SQRT(ELAM)*BETA)      

      IF (IEND.EQ.0) THEN
         IF (ABS(LAMBDM).GT.1.D0) THEN
             LAMBDA = LAMBDP 
         ELSE
             LAMBDA = LAMBDM
         END IF

         PDY = Y22B-LAMBDA   

      ELSE
         IF (ABS(LAMBDM).GT.1.D0) THEN
             LAMBDA = LAMBDM 
         ELSE
             LAMBDA = LAMBDP
         END IF

         PDY = Y22B-LAMBDA

      END IF

      IF (ABS(Y)+ABS(PDY).GT.1.D-3) RETURN

      Y = LAMBDA+SIN(SQRT(ELAM-ALPHA)*BETA)
     &  *SIN(SQRT(ELAM)*BETA)/SQRT(ELAM-ALPHA)
     & - COS(SQRT(ELAM)*BETA)*COS(SQRT(ELAM-ALPHA)*BETA)

      PDY = -SQRT(ELAM)*COS(SQRT(ELAM-ALPHA)*BETA)
     &    *SIN(SQRT(ELAM)*BETA) - SQRT(ELAM-ALPHA)
     &    *SIN(SQRT(ELAM-ALPHA)*BETA)*COS(SQRT(ELAM)*BETA)

      WRITE (11,FMT=*) 'Y,PDY:',Y,PDY


      RETURN

120   IF (IEND.EQ.0) THEN
        ALF = ELAM+CMPLX(0.5D0,0.D0)
        IF (DBLE(ALF).GE.0.D0) THEN
            ALF = CMPLX(0.5D0,0.D0)+ALF
        ELSE
            ALF = CMPLX(0.5D0,0.D0)-ALF
        END IF
        Z = XI(C12,CMPLX(COS(NU),SIN(NU))*X)
        QL0 = C12**2 - 1.D0/(1.D0+Z)+0.25D0/((C12*(1.D0+Z))**2)
        DQL0 = 1.D0/(1.D0+Z)**2 - 0.5D0/((C12**2)*((1.D0+Z)**3))
C        WRITE (9,FMT=*) 'EPS,NU (L):',X,NU
        A1 = 0.5D0/ALF
        A2 = (A1-C12**2)/(4.D0*ALF+2.D0)
        A3 = (A2-A1*C12**2)/(6.D0*(ALF+1.D0))
        A4 = (A3-A2*C12**2)/(8.D0*ALF+12.D0)
        Y = 1.D0 + Z*(A1 + Z*(A2 + Z*(A3 + A4*Z))) 
        PDY = (CMPLX(COS(NU),SIN(NU))/Z)*(ALF + Z*((ALF+1.D0)*A1+ 
     &          Z*((ALF+2.D0)*A2 + Z*((ALF+3.D0)*A3 +
     &          Z*(ALF+4.D0)*A4))))
        PDY = (PDY/Y + 0.25D0*DQL0/QL0)/SQRT(QL0)
        Y = CMPLX(1.D0,0.D0)
      ELSE
C        Z = XI(C12,X*CMPLX(COS(NU),SIN(NU)))
C        QL0 = C12**2 - 1.D0/(1.D0+Z)+0.25D0/((C12*(1.D0+Z))**2)
C        DQL0 = 1.D0/(1.D0+Z)**2 - 0.5D0/((C12**2)*((1.D0+Z)**3))
C        QLOC = C12*C12-1.D0/Z-(ELAM*ELAM+ELAM)/(Z*Z)
C        DQLOC = 1.D0/(Z*Z)+2.D0*(ELAM*ELAM+ELAM)/(Z*Z*Z)
C        SQLOC = SQRT(QLOC)
C        WRITE (9,FMT=*) 'SQLOC:',SQLOC
C        WRITE (9,FMT=*) 'EPS,NU (R):',EPS,NU
C        IF (DBLE(SQLOC)*DBLE(C12).LT.0.D0) SQLOC = -SQLOC
C        Y = CMPLX(1.D0,0.D0)
C        PDY = CMPLX(COS(NU),SIN(NU))*(CMPLX(0.D0,1.D0)*SQLOC
C     &        -0.25D0*DQLOC/QLOC)
      END IF
      PDY = CMPLX(0.D0,1.D0)
      Y = CMPLX(1.D0,0.D0)
      RETURN

130   IF (IEND.EQ.0) THEN
        Y = CMPLX(0.D0,0.D0)
        PDY = CMPLX(1.D0,0.D0)
      ELSE
        Y = CMPLX(1.D0,0.D0)
        PDY = SQRT(-ELAM)
        IF (DBLE(PDY).GT.0.D0) PDY = -PDY
      END IF
      RETURN

140   CONTINUE
      Z = EPS+CMPLX(COS(NU),SIN(NU))*X
      IF (IEND.EQ.0) THEN
        IF (PE.GT.2) THEN
           Y = CMPLX(0.D0,0.D0)
           PDY = CMPLX(1.D0,0.D0)
        ELSE IF (PE.EQ.2) THEN
           Y = CMPLX(1.D0,0.D0)
           PDY = 0.5D0*SQRT(5.D0+4.D0*ELAM*(ELAM+1.D0))
           IF (DBLE(PDY).LT.0.D0) PDY = -PDY
           PDY = PDY+0.5D0
        ELSE
           Y = CMPLX(1.D0,0.D0)
           PDY = ELAM+CMPLX(1.D0,0.D0)
        END IF
      ELSE
        Y = CMPLX(1.D0,0.D0)
        PDY = CMPLX(0.D0,1.D0)*C14*
     &        CMPLX(COS(NU),SIN(NU))*(1.D0 
     &        - 0.5D0/((C14**2)*(Z**PE))
     &        -0.5D0*(ELAM*ELAM+ELAM)/(C14*C14*Z*Z))
C        WRITE (7,FMT=*) 'Kappa,NU:',C14,NU
      END IF
      RETURN

150   CONTINUE
      IF (IEND.EQ.0) THEN
        Y = CMPLX(1.D0,0.D0)
        PDY = ELAM + CMPLX(0.5D0,0.D0)
        IF (DBLE(PDY).LT.0.D0) PDY = -PDY
      ELSE
        Y = CMPLX(0.D0,0.D0)
        PDY = CMPLX(1.D0,0.D0)
      END IF

160   CONTINUE
      IF (IEND.EQ.0) THEN
        Y = CMPLX(1.D0,0.D0)
        IF (PE.LT.2) THEN
            PDY = ELAM + CMPLX(0.5D0,0.D0)
            IF (DBLE(PDY).LT.0.D0) PDY = -PDY
        ELSE IF (PE.EQ.2) THEN
            PDY = SQRT((ELAM+0.25D0)**2+C161)
            IF (DBLE(PDY).LT.0.D0) PDY = -PDY
        ELSE
             Y = CMPLX(0.D0,0.D0)
             PDY = CMPLX(1.D0,0.D0)
        END IF
      ELSE
        Y = CMPLX(0.D0,0.D0)
        PDY = CMPLX(1.D0,0.D0)
      END IF

      RETURN

170   CONTINUE
      IM = CMPLX(0.D0,1.D0)
      IF (IEND.EQ.0) THEN
            Y = CMPLX(0.D0,0.D0)
            PDY = CMPLX(1.D0,1.D0)
            KAPPA = C162
            CEE = C161
            IX = CMPLX(X,0.D0)
            CALL LEFTBC(IX,KAPPA,CEE,ELAM,Y,PDY)
      ELSE

           Y = CMPLX(1.D0,0.D0)
           IX = CMPLX(X,0.D0)
           B1 = -(C161/(2.D0*C162*IM) 
     &            + C161/(4.D0*(C162**2))-(ELAM*ELAM+ELAM))/
     &                (2.D0*C162*IM)
           B2 = (2.D0-C161/(C162*IM))*B1/
     &           (4.D0*C162*IM)
           B3 = (B1*B1+6.D0*B1-C161*2.D0*B1/(C162*IM))/
     &           (6.D0*C162*IM)
           PDY = IM*C162+C161/(2.D0*C162*IM*IX)-B1/(IX**2)
     &           -2.D0*B2/(IX**3)-3.D0*B3/(IX**4)
           
      END IF
      RETURN

180   CONTINUE
      IF (IEND.EQ.0) THEN
            Y = CMPLX(0.D0,0.D0)
            PDY = CMPLX(1.D0,1.D0)
      ELSE
            Y = CMPLX(1.D0,1.D0)
            PDY = CMPLX(0.D0,0.D0)
      END IF
      RETURN

190   CONTINUE
      IF (IEND.EQ.0) THEN

            Q = C161/(X**4)-C162**2+(ELAM+ELAM**2)/(X**2)
            DQ = -4.D0*C161/(X**5) -2.D0*(ELAM+ELAM**2)/(X**3)
            D2Q = -20.D0*C161/(X**6) - 6.D0*(ELAM+ELAM**2)/(X**4)

            DS0 = SQRT(Q)
            DS1 = -DQ/(4.D0*Q)           
            DS2 = (D2Q/(8.D0*Q)-5.D0*(DQ**2)/(32.D0*Q*Q))/DS0
 
            Y = CMPLX(1.D0,0.D0)
            PDY = DS0 + DS1 + DS2
C            PDY = SQRT(C161/(X**4))

      ELSE

           Y = CMPLX(1.D0,0.D0)
           IX = CMPLX(X,0.D0)
           PDY = IM*C162-IM*(ELAM+ELAM**2)/(2.D0*C162*(IX**2))

      END IF

      RETURN

200   IF (IEND.EQ.0) THEN
         R = CMPLX(X,0.D0)
         CALL FERMBC(R,ELAM,C163,C164,C161,C162,Y,PDY)
      ELSE
         Y = CMPLX(1.D0,0.D0)
         IX = CMPLX(X,0.D0)
         IM = CMPLX(0.D0,1.D0)
         PDY = IM*C162+C161/(2.D0*C162*IM*IX)

           B1 = (ELAM*ELAM+ELAM)/(2.D0*C162*IM)
           B2 = 2.D0*B1/(4.D0*C162*IM)
           B3 = (B1*B1+6.D0*B1)/(6.D0*C162*IM)
           PDY = IM*C162-B1/(IX**2)
     &           -2.D0*B2/(IX**3)-3.D0*B3/(IX**4)

      END IF
C      WRITE (19,FMT=*) 'IEND,ELAM,Y,PDY',IEND,ELAM,Y,PDY 

      RETURN

210   CALL GETRT(X,C161,C162,C163,C164,ELAM,R,DR,D2R)
      IF (IEND.EQ.0) THEN
         CALL FERMBC(R,ELAM,C163,C164,C161,C162,Y,PDY)
         PDY = DR*PDY
      ELSE
         Y = CMPLX(1.D0,0.D0)
         IX = R
         IM = CMPLX(0.D0,1.D0)
         PDY = IM*C162+C161/(2.D0*C162*IM*IX)

         B1 = (ELAM*ELAM+ELAM)/(2.D0*C162*IM)
         B2 = 2.D0*B1/(4.D0*C162*IM)
         B3 = (B1*B1+6.D0*B1)/(6.D0*C162*IM)
         PDY = IM*C162-B1/(IX**2)
     &           -2.D0*B2/(IX**3)-3.D0*B3/(IX**4)
         PDY = DR*PDY
      END IF 
      RETURN

220   IF (IEND.EQ.0) THEN
         Y = CMPLX(1.D0,0.D0)
         PDY = (ELAM+1.D0)*Y/X
      ELSE
         Y = CMPLX(1.D0,0.D0)
         PDY = (ELAM+1.D0)*Y/X
      END IF
C      WRITE (19,FMT=*) 'IEND,ELAM,Y,PDY',IEND,ELAM,Y,PDY 
      RETURN

      END


C -------------------------------------------------------------------

      DOUBLE COMPLEX FUNCTION MYSQRT(Z)
      DOUBLE COMPLEX Z
      DOUBLE PRECISION X,Y,THETA
C Computes the square root whose argument lies in [Pi,2.Pi).
      X = DBLE(Z)
      Y = AIMAG(Z)
C Get the argument of Z:
      THETA = ATAN2(Y,X)
      IF (THETA.LT.0.D0) THEN 
C Add 2.Pi to the argument of Z 
          THETA = THETA + 8.D0*ATAN(1.D0)
      END IF
C Arg(Z) now lies in [0,2.Pi)
      THETA = THETA/2.D0
      MYSQRT = -((X*X+Y*Y)**0.25D0)*CMPLX(COS(THETA),SIN(THETA))
      RETURN
      END

C ------------------------------------------------------------------

      SUBROUTINE SLSMPC(X,ELAM,P,Q1,Q0)

C     .. Scalar Arguments ..
      DOUBLE COMPLEX ELAM,Q0,Q1
      DOUBLE PRECISION P,X
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC CMPLX,AIMAG
C     ..
C
C The comparison problem is -y'' + (1+Re(lambda))y = \lambda y
C 
C NB Re(lambda) is constant on a vertical line in the complex plane.
C 
      P = 1.D0
      Q1 = CMPLX(0.D0,0.D0)
      Q0 = CMPLX(1.D0,-AIMAG(ELAM))

      RETURN
      END

C ------------------------------------------------------------------


      SUBROUTINE SLSMBC(ISING,X,ELAM,Y,PDY)

C     .. Scalar Arguments ..
      DOUBLE COMPLEX ELAM,Y,PDY
      DOUBLE PRECISION X
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC CMPLX
C     ..

C
C Only Dirichlet BCs are currently implemented for the comparison problem
C
      Y = CMPLX(0.D0,0.D0)
      PDY = CMPLX(1.D0,0.D0)

      RETURN
      END

C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      SUBROUTINE BCXTRP(C,K,XO,XBIG,TOL,ELL,U,DU)
      DOUBLE COMPLEX  C,K,ELL,U,DU,U1,DU1,U2,DU2,U3,DU3,
     &            U4,DU4
      DOUBLE PRECISION XO,XBIG,TOL

      CALL BCCALC(C,K,XO,XBIG,TOL,ELL,U,DU)
      CALL BCCALC(C,K,XO,2.D0*XBIG,TOL,ELL,U1,DU1)
      CALL BCCALC(C,K,XO,4.D0*XBIG,TOL,ELL,U2,DU2)

      U3 = 2.D0*U1-U
      DU3 = 2.D0*DU1-DU
      U4 = 2.D0*U2-U1
      DU4 = 2.D0*DU2-DU1

      U = (4.D0*U4-U3)/3.D0
      DU = (4.D0*DU4-DU3)/3.D0

      RETURN
      END


      SUBROUTINE BCCALC(C,K,XO,XBIG,TOL,ELL,U,DU)
      IMPLICIT NONE
      DOUBLE COMPLEX C,K,ELL,U,DU,U3,DU3
      DOUBLE PRECISION XO,XBIG,TOL
      DOUBLE COMPLEX CO,KO,ELLO
      COMMON /PARAMS/CO,KO,ELLO
      INTEGER NEQ,LENWRK
      PARAMETER(NEQ=4,LENWRK=32*NEQ)
      DOUBLE PRECISION YGOT(NEQ),YPGOT(NEQ),YMAX(NEQ),THRES(NEQ),
     &                 WORK(LENWRK)

      DOUBLE PRECISION TWANT,TGOT
      INTEGER I,IFAIL,METHOD

      EXTERNAL F,X02AMF
      DOUBLE PRECISION X02AMF

      CO = C
      KO = K
      ELLO = ELL

      U3 = CMPLX(0.D0,0.D0)
      DU3 = CMPLX(0.D0,0.D0)

      YGOT(1) = DBLE(U3)
      YGOT(2) = AIMAG(U3)
      YGOT(3) = DBLE(DU3)
      YGOT(4) = AIMAG(DU3)
      TWANT = XO
      TGOT = XBIG
      IFAIL = 0
      DO 10 I = 1,NEQ
         THRES(I) = SQRT(X02AMF(0.D0))
10    CONTINUE
      METHOD = 3 

      CALL D02PVF(NEQ,XBIG,YGOT,XO,TOL,THRES,METHOD,'U',.FALSE.,
     &            0.D0,WORK,LENWRK,IFAIL)

      IFAIL = 1

20    CALL D02PCF(F,TWANT,TGOT,YGOT,YPGOT,YMAX,WORK,IFAIL)
      IF (IFAIL.EQ.3) GOTO 20

      U3 = CMPLX(YGOT(1),YGOT(2))
      DU3 = CMPLX(YGOT(3),YGOT(4))

      U = U3 + CMPLX(0.D0,1.D0)*K*(XO-0.5D0*C*LOG(XO)/(K**2))
      DU = DU3 + CMPLX(0.D0,1.D0)*K*(1.D0-0.5D0*C/(XO*(K**2)))

      RETURN
      END

      SUBROUTINE F(X,Y,YP)
      IMPLICIT NONE
      DOUBLE COMPLEX C,K,ELL
      COMMON /PARAMS/C,K,ELL      
      DOUBLE PRECISION X
      DOUBLE PRECISION Y(4),YP(4)
      DOUBLE COMPLEX U3,DU3,G(2),I

      U3 = CMPLX(Y(1),Y(2))
      DU3 = CMPLX(Y(3),Y(4))
      I = CMPLX(0.D0,1.D0)

      G(1) = DU3
      G(2) = -I*C/(2.D0*K*X*X)+C*C/(4.D0*K*K*X*X)
     &       +ELL*(ELL+1.D0)/(X*X)-DU3**2
     &       -2.D0*DU3*I*K*(1.D0-C/(2.D0*K*K*X))

      YP(1) = DBLE(G(1))
      YP(2) = AIMAG(G(1))
      YP(3) = DBLE(G(2))
      YP(4) = AIMAG(G(2))
C      WRITE (6,FMT=*) 'YP',YP

      RETURN
      END

C -------------------------------------------------------------------

      SUBROUTINE LEFTBC(X,K,C,ELAM,Y,PDY)
      IMPLICIT NONE
      DOUBLE COMPLEX K,C,ELAM,Y,PDY
      DOUBLE COMPLEX X
      INTEGER N
      PARAMETER (N=40)
      DOUBLE COMPLEX A(0:N),ALFA
      INTEGER I

      ALFA = ELAM+CMPLX(1.D0,0.D0)

      A(0) = CMPLX(1.D0,0.D0)
      A(1) = C/(2.D0*ALFA)
      DO 10 I = 2,N
        A(I) = (C*A(I-1)-K*K*A(I-2))/(2.D0*CMPLX(I,0.D0)*ALFA
     &          +CMPLX(I*(I-1.D0),0.D0))
10    CONTINUE

      Y = A(N)
      DO 20 I = N,1,-1
         Y = Y*X+A(I-1)
20    CONTINUE
C      Y = Y*(X**ALFA)

      PDY = A(N)*(N+ALFA)
      DO 21 I = N,1,-1
        PDY = PDY*X+A(I-1)*(I-1+ALFA)
21    CONTINUE
C      PDY = PDY*(X**(ALFA-1))
      PDY = PDY/X

      WRITE (17,FMT=*) 'ELAM,Y,PDY:',ELAM,Y,PDY
      WRITE (17,FMT=*) 'A(0),A(N)X**N:',A(0),A(N)*(X**N)
      WRITE (17,FMT=*) 'C,KAPPA:',C,K

      RETURN
      END

C -----------------------------------------------------------------

      SUBROUTINE FERMBC(X,ELAM,A,B,C,KAPPA,Y,PDY)
      IMPLICIT NONE
      INTEGER I,N,M
      PARAMETER (N=40)
      DOUBLE COMPLEX X  
      DOUBLE COMPLEX ELAM,A,B,C,KAPPA,AL,BL
      DOUBLE COMPLEX AN(0:N),S(1:N),BEE(0:N)
      DOUBLE COMPLEX ALFA
      DOUBLE COMPLEX Y,PDY

C Correction to silly mistake made while writing routine:
      AL = -A
      BL = -B
C

      ALFA = ELAM + CMPLX(1.D0,0.D0) 
C      ALFA = -ELAM

      AN(0) = CMPLX(1.D0,0.D0)
      AN(1) = 0.5D0*C*AN(0)/ALFA


      S(1) = AN(0)
      BEE(1) = AN(0)
      S(2) = AN(1) + AL*AN(0)
      BEE(2) = AN(1) + AL*AN(0)

      AN(2) = (C*BEE(2)-(KAPPA**2)*AN(0))/(4.D0+2.D0*(2.D0*ALFA-1.D0))

      DO 10 M = 3,N
     
         S(M) = AL*S(M-1)+AN(M-1)
         BEE(M) = BL*BEE(M-2)+S(M)
         AN(M) = (C*BEE(M)-(KAPPA**2)*AN(M-2))/(M*M+M*(2.D0*ALFA-1.D0))

C      WRITE (17,FMT=*) 'M,AN(M),BEE(M),S(M):',M,AN(M),BEE(M),S(M)

10    CONTINUE
  


      Y = AN(N)
      PDY = (N+ALFA)*AN(N)
      DO 20 I = N,1,-1
         Y = Y*X + AN(I-1)
         PDY = PDY*X + (I-1+ALFA)*AN(I-1)
20    CONTINUE
      WRITE (17,FMT=*) 'N*A(N)*X**N:',N*AN(N)*(X**N)
      PDY = PDY/X

      RETURN
      END
