\begin{verbatim}
c
c********************************************************************
c
c Eigenvalues for scalar sturm-liouville equations with arbitrary
c separated or coupled boundary conditions.
c
c********************************************************************
c This is a program based on the theory of Atkinson and the
c work of Atkinson, Krull, Lief, Zettl which will compute the
c eigenvalues of matrix equations.
c
c The equation we are working with is:
c
c         -(py')' + qy = lambda*wy
c
c       C11*y(a) + C12*py'(a) + D11*y(b) + D12*py'(b) = 0
c       C21*y(a) + C22*py'(a) + D21*y(b) + D22*py'(b) = 0
c
c written by H. Dwyer  Oct. 91
c modified Jun. 95
c modified Aug. 95
c--------------------------------------------------------------------
c VARIABLE DECLARATIONS
c
      double precision eigtol,lo,midpt,hi,span,a,b,cntr,pi,lobnd
c      
      double complex b1(2,2),b2(2,2),cmplxi,
     * cmplx0,cmplx1
      double complex capba(2,2),capta(2,2),bfctr(2,2),tmp1(2,2),
     * tmp2(2,2),vcapba(2)
c      
      integer neig,targt,nsplit,locnt,midcnt,hicnt
      integer track,qda(2),quad,offset
c
c--------------------------------------------------------------------
c Description of Variables
c
c b1,b2       : complex constant 2x2 matrices which define the 
c               boundary cond. for the 2x2 prob
c eigtol      : the user supplied tolerance for the result
c lo,midpt,hi : low end, midpoint, and high end of the bisection 
c               bracketing interval
c span        : the length of the initial bracketing interval, 
c               used to compute the number of bisections required
c a,b,cntr    : lower,upper endpts, midpt of the problem interval
c pi          : the irrational constant, about 3.14159
c cmplx0      : the complex constant 0+0i
c cmplx1      : the complex constant 1+0i
c cmplxi      : the complex constant 0+1i
c lobnd       : a lower bound on the spectrum
c capba,capta : initial values of the matrices Theta and B at "cntr"
c vcapba      : the eigenvectors of B(cntr)
c bfctr       : the constant conversion factor to get B from Theta
c tmp1,tmp2   : temporary matrix work space
c neig        : the user supplied number of the desired eigenvalue
c offset      : offset used to find the target index
c targt       : the target index jump corresponding to eigenvalue 
c               neig
c nsplit      : the number of bisections required to get an 
c               interval of length < tol
c locnt,hicnt : the index count for the low and high bracketing 
c               values
c midcnt      : the index count for the midpoint of the bisection 
c               bracket
c bndcnt      : the index count for the spectrum lower bound, lobnd
c track       : an integer flag to enable/disable diagnostic output
c qda         : the quadrants of the eigenvalues of B(cntr)
c quad        : integer quadrant function
c
c--------------------------------------------------------------------
c COMMON BLOCK ASSIGNMENTS
c
      common/blk1/pi,a,b,cntr,b1,b2,track
      common/blk2/hi,lo,lobnd,targt,locnt,hicnt
      common/blk3/capta,bfctr,qda,vcapba
c
c--------------------------------------------------------------------
c Data initializations
c
      pi=4.0d0*datan(1.0d0)
      cmplx0=dcmplx(0.0d0,0.0d0)
      cmplx1=dcmplx(1.0d0,0.0d0)
      cmplxi=dcmplx(0.0d0,1.0d0)
c
c--------------------------------------------------------------------
c User Input
c
c --- We ask the user to specify which eigenvalue they want to 
c     compute and the tolerance that they want it computed to
c
      write(6,*)' '
      write(6,*)'Matrix eigenvalue solver for regular scalar'
      write(6,*)'problems using simple 4th order Runge-Kutta Method'
      write(6,*)' '
      write(6,*)'Eigenvalues are numbered beginning with #1'
      write(6,*)'That is, the first eigenvalue is #1, not #0'
      write(6,*)' '
      write(6,*)'Which eigenvalue do you wish to compute ?'
      read(5,*) neig
   50 write(6,*)'What tolerance would you like it to ?'
      read(5,*) eigtol
      if(eigtol.le.0.0d0) goto 50
c
c --- begin by setting up the interval and initializing parameters
c
c --- check to see if the user wants to watch the bisection and 
c     other information printouts
   51 write(6,*)'Do you wish to see more informative output ?'
      write(6,*)'    Enter 1 for informative output'
      write(6,*)'    enter 0 to keep it brief'
      read(5,*)track
      if((track.ne.0).and.(track.ne.1)) goto 51
c
c--------------------------------------------------------------------
c --- call the user supplied routine which sets up the problem
c
      call setup 
c
c--------------------------------------------------------------------
c --- we compute the initial values of the capt matrix at cntr, and a
c     conversion factor to get capb from capt
c     the "boundary conditions" at the center are always the same,
c     we can compute capt once (by hand) and simply assign the result
c
      capta(1,1)=cmplx0
      capta(1,2)=-cmplx1
      capta(2,1)=-cmplx1
      capta(2,2)=cmplx0
c
c                                      -1
c --- we will compute   bfctr = (B1-iB2) *(B1+iB2)
c
c --- we form tmp2=B1-iB2
c
      do 120 i=1,2
        do 120 j=1,2
          tmp2(i,j)=b1(i,j)-cmplxi*b2(i,j)
  120 continue
c                                        -1
c --- call routine to compute tmp1 = tmp2
c
      call invert(tmp2,tmp1,ifail)
c
c --- check to be sure that the INVERT routine succeeded
c     if not, go to the error handling part of our program
c
      if(ifail.ne.0) goto 1100
c
c --- we form tmp2=B1+iB2
c
      do 130 i=1,2
        do 130 j=1,2
          tmp2(i,j)=b1(i,j)+cmplxi*b2(i,j)
  130 continue
c                       
c --- we form bfctr=tmp1 *tmp2
c
      call matmul(tmp1,0,tmp2,0,bfctr)
c
c --- we form B(cntr) = bfctr * T(cntr)
c     we use a routine to get the eigenvalues of B(cntr)
c     then we get the quadrants for the eigenvalues of B(cntr)
c
      call matmul(bfctr,0,capta,0,capba)      
      call eigval(capba,vcapba(1),vcapba(2))
c
      do 140 i=1,2
        qda(i)=quad(vcapba(i))
  140 continue
c
c--------------------------------------------------------------------
c --- select the index count that we wish to match - that is, the
c     number of times an eigenvalue of capb has passed through 1
c     with an offset calculated so that the labels are correct
c
      call getoff(offset)
c
c --- compute the target index value using the offset and the
c     user supplied neig
c
      targt=neig+offset
c
c--------------------------------------------------------------------
c --- find a pair of bracketing values for the eigenvalue sought
c
      call brkt
c
      if(track.eq.1) then
        write(6,*)' '
        write(6,*)'The target index is ',targt 
        write(6,*)' '
        write(6,*)'The lo value is ',lo
        write(6,*)'The hi value is ',hi
        write(6,*)'The lo index is ',locnt
        write(6,*)'The hi index is ',hicnt
      endif
c
c--------------------------------------------------------------------
c --- refine our bracketing values using a simple bisection
c     algorithm until the difference between the hi and lo values 
c     is less than the tolerance given by the user.
c
c --- compute the number of bisection steps required to refine the
c     values as required
c
      span=hi-lo
      nsplit=int(log(span/eigtol)/log(2.0)) + 1
c
      write(6,*)' '
      write(6,*)'Beginning bracket refinement'
      if(track.eq.1) then
        write(6,*)' '
        write(6,*)'We will refine the bracket by simple bisection'
        write(6,*)'until the eigenvalue is determined to within the'
        write(6,*)'tolerance you have specified.'
        write(6,*)' '
        write(6,*)'The number of bisections required is',nsplit
      endif
c
c --- we now perform the bisections
c
      do 200 indx1=1,nsplit
        write(6,*)'    Bisections remaining:',(nsplit-indx1)
c
c ---   compute the midpoint and get the index count for that value
c
        midpt=(hi+lo)/2.0d0
        midcnt=kount(midpt)
c
c
c ---   if the index count for the midpoint is less then the target
c       count then we have a new lo, otherwise we have a new hi
c
        if(midcnt.lt.targt) then
          lo=midpt
          locnt=midcnt
        else
          hi=midpt
          hicnt=midcnt
        endif
c
        if(track.eq.1) then
          write(6,*)' Lower endpoint ',lo,' has an index of ',locnt
          write(6,*)' Upper endpoint ',hi,' has an index of ',hicnt
          write(6,*)' '
        endif
c
  200 continue
c
c--------------------------------------------------------------------
c --- we are done with the bisections, so we compute the midpoint
c     one last time, and report the result
c
      midpt=(lo+hi)/2.0d0
      write(6,*)' %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% '
      write(6,901) neig
      write(6,902) midpt
      write(6,903) lo,hi
      write(6,904) (hicnt-locnt)
  901 format(' Eigenvalue number ',i3)
  902 format(' Best guess for the value is ', f12.6)
  903 format(' Interval est. (',f12.6,' ,',f12.6,' )')
  904 format(' Multiplicity is ',i4)
      goto 9999
c
c--------------------------------------------------------------------
c    Error handling section for main routine
c
 1100 continue
c
c --- the routine for finding inverses of complex matrices has 
c     failed
      write(6,*)' The routine INVERT has failed'
      write(6,*)' Check B.C. matrices for proper values.'
      goto 9999
c
c
 9999 continue
      end                                                          
c
c********************************************************************
c                                    
c  CHECKBC  
c
c********************************************************************
c
c This routine checks that the boundary conditions supplied by the
c user satisfy the requirements for a self-adjoint problem
c 
      subroutine checkbc(bc1,bc2,flag)
c
      double complex bc1(2,2),bc2(2,2),tmp1,tmp2,cmplx0
      real tmp3,tmp4
      integer flag
c
      cmplx0=dcmplx(0.0d0,0.0d0)
c
c --- we begin by assuming that we will succeed
c
      flag=0
c 
      write(6,*)' '
      write(6,*)'Checking Boundary conditions'
      write(6,*)' '
c                         ___   ___           ___   ___
c --- requirement #1  C11*C12 - C11*C12 = D11*D12 - D11*D12
c
      tmp1=bc1(1,1)*dconjg(bc1(1,2))-dconjg(bc1(1,1))*bc1(1,2)
      tmp2=bc2(1,1)*dconjg(bc2(1,2))-dconjg(bc2(1,1))*bc2(1,2)
      tmp3= zabs(tmp1-tmp2)
      if((tmp3+1.0).ne.1.0) flag=1
c                         ___   ___           ___   ___
c --- requirement #2  C21*C22 - C21*C22 = D21*D22 - D21*D22
c
      tmp1=bc1(2,1)*dconjg(bc1(2,2))-dconjg(bc1(2,1))*bc1(2,2)
      tmp2=bc2(2,1)*dconjg(bc2(2,2))-dconjg(bc2(2,1))*bc2(2,2)
      tmp3= zabs(tmp1-tmp2)
      if((tmp3+1.0).ne.1.0) flag=2
c                         ___   ___           ___   ___
c --- requirement #3  C21*C12 - C11*C22 = D21*D12 - D11*D22
c
      tmp1=bc1(2,1)*dconjg(bc1(1,2))-dconjg(bc1(1,1))*bc1(2,2)
      tmp2=bc2(2,1)*dconjg(bc2(1,2))-dconjg(bc2(1,1))*bc2(2,2)
      tmp3= zabs(tmp1-tmp2)
      if((tmp3+1.0).ne.1.0) flag=3
c
c --- requirement #4  rank(C|D) = 2
c
      tmp3=zabs(bc1(1,1))+zabs(bc1(1,2))+zabs(bc2(1,1))+
     2     zabs(bc2(1,2))
      tmp4=zabs(bc1(2,1))+zabs(bc1(2,2))+zabs(bc2(2,1))+
     2     zabs(bc2(2,2))
      if(((tmp3+1.0).eq.1.0).or.((tmp4+1.0).eq.1.0)) then
        flag=4
      else
        if(bc1(1,1).ne.cmplx0) then
          tmp1=bc1(2,1)/bc1(1,1)
        elseif(bc1(1,2).ne.cmplx0) then
          tmp1=bc1(2,2)/bc1(1,2)
        elseif(bc2(1,1).ne.cmplx0) then
          tmp1=bc2(2,1)/bc2(1,1)
        else
          tmp1=bc2(2,2)/bc2(1,2)
        endif
        tmp3=0.0
        tmp3=tmp3+zabs(bc1(2,1)-tmp1*bc1(1,1))
        tmp3=tmp3+zabs(bc1(2,2)-tmp1*bc1(1,2))
        tmp3=tmp3+zabs(bc2(2,1)-tmp1*bc2(1,1))
        tmp3=tmp3+zabs(bc2(2,2)-tmp1*bc2(1,2))
        if((tmp3+1.0).eq.1.0) flag=4
      endif
      return
      end
c
c
c********************************************************************
c                                    
c  BRKT  
c
c********************************************************************
c
c this subroutine will find a bracketing pair of values
c which enclose the desired eigenvalue
c
c the user may supply a guess if they wish, otherwise
c we start with lobnd
c
      subroutine brkt
c
      double complex b1(2,2),b2(2,2)
      double precision hi,lo,lobnd,test,step,pi,a,b,cntr
      integer targt,locnt,hicnt,maxi,track
      logical lofnd,hifnd
c
      common/blk1/pi,a,b,cntr,b1,b2,track
      common/blk2/hi,lo,lobnd,targt,locnt,hicnt
c
c --- initialize the logical variables which tell 
c     whether we have a bracket
c
      lofnd=.false.
      hifnd=.false.
c
c --- initialize the step size for our search, and the
c     max number of attempts that are allowed
c
      step=4.0d0
      maxi=20
      if(track.eq.1) then
        write(6,*)' '
        write(6,*)'We begin by looking for two values which bracket'
        write(6,*)'the desired eigenvalue.  This is done using a'
        write(6,*)'simple step-by-step search, with a stepsize'
        write(6,*)'of:',step
        write(6,*)' '
        write(6,*)'The maximum number of steps which will be allowed'
        write(6,*)'before the program "gives up" is:',maxi
        write(6,*)' '
        write(6,*)'A good guess, which will be used as a starting'
        write(6,*)'point for our search, would be helpful.'
        write(6,*)' '
      endif
c
c --- we begin by asking the user if they have a guess
c
   10 write(6,*)'Do you have a guess for the eigenvalue ?'
      write(6,*)'     0=no, 1=yes'
      read(5,*)ians
      if((ians.gt.1).or.(ians.lt.0))goto 10
c
c --- if they have a guess, go get it
c     if not, use lobnd
c
      if(ians.eq.0) then
        test=lobnd
      else
        write(6,*)'Input guess value.'
        read(5,*) test
      endif
c
c --- the search begins, and continues until both ends are
c     found, or until we have tried too many times
c
      write(6,*)' '
      write(6,*)'Beginning search for bracketing values'
      do 100 i=1,maxi
        write(6,*)'    Bracket search step ',i
        itest=kount(test)
        if(itest.lt.targt) then
c
c         test is a lower bracketing value
c
          lo=test
          if(track.eq.1) 
     *      write(6,*)'    ',lo,' is a lower bracketing value'
          locnt=itest
          lofnd=.true.
          test=test+step
c
        else
c
c         test is a high bracketing value
c
          hi=test
          if(track.eq.1) 
     *      write(6,*)'    ',hi,' is a upper bracketing value'
          hicnt=itest
          hifnd=.true.
          test=test-step
c
        endif
c
c ---   check to see if we have both ends
c       if not, go back and try again
c
        if(lofnd.and.hifnd) goto 900
  100 continue
c
c --- if we get here, we did not find a bracket in the max number
c     of attempts
      write(6,*)' After ',maxi,' attempts, bracket search'
      if(.not.hifnd) 
     2  write(6,*)' Failed to find an upper bracket endpt.'
      if(.not.lofnd) 
     2  write(6,*)' Failed to find a lower bracket endpt.'
      write(6,*)' Last test point was ',test
      stop
c
c     we are done
c
  900 continue
      return
      end
c
c
c********************************************************************
c
c  GETOFF
c
c********************************************************************
c
c this is a subroutine which will compute the offset 
c required so that the eigenvalues are labelled
c correctly
c
      subroutine getoff(offset)
c
      double precision lobnd,lo,hi
      integer offset,targt,locnt,hicnt
c
      common/blk2/hi,lo,lobnd,targt,locnt,hicnt
c
  10  write(6,*)'Is the spectrum for this problem bounded below?'
      write(6,*)'    0=no, 1=yes'
      read(5,*)ians
      if((ians.gt.1).or.(ians.lt.0)) goto 10
c
c --- once we know which case we are in, we know what to do
c
      if(ians.eq.0) then
c 
c       in the unbounded case, we compute the count for
c       a "lobnd" of zero, and use it to get the offset
c
        lobnd=0.0d0
        offset=kount(lobnd)+1
c
      else
c
c       in the bounded case, we ask the user for a lower bound,
c       and use that to find the offset
c
        write(6,*)'Please input a lower bound.'
        read(5,*)lobnd
        offset=kount(lobnd)
c    
      endif
      return
      end
c
c
c********************************************************************
c
c  KOUNT
c
c this function finds the index count for a specified lambda
c
c********************************************************************
c
      function kount(lamda)
c
      integer kount,mesh,icnt,track,qda(2),qdold(2),
     * qd(2),quad,submsh
      double precision lamda,pi,a,b,h,x,x1,x2
      double complex capta(2,2),bfctr(2,2),
     * capt(2,2),capb(2,2),vcapba(2),
     * h12,h16,h13
      double complex cmplxh,cmplx1,vcapb(2),b1(2,2),b2(2,2)
      double complex w1(2,2),w2(2,2),w3(2,2),
     * t1(2,2),t2(2,2),t3(2,2),t4(2,2)
c
      common/blk1/pi,a,b,cntr,b1,b2,track
      common/blk3/capta,bfctr,qda,vcapba
c
c--------------------------------------------------------------------
c Description of Variables
c
c mesh           : number of locations in [cntr,b] where eigenvalues 
c                  are computed, and tracked as they move around 
c                  the circle
c submsh         : number of nodes between adjacent mesh pts, where
c                  integration is done, but no eigenvalues computed 
c h              : stepsize between submesh points
c qd, qda, qdold : quadrants of eigenvalues for current B matrix, for
c                B matrix at endpoint "cntr", and at the prev mesh pt
c lamda          : the real value whose index count we want (yes, I
c                  know it's spelled wrong)
c pi             : the irrational constant, about 3.1416
c a,b            : the endpoints of the interval
c b1,b2          : complex matrices defining boundary conditions
c track          : integer flag to enable/disable diagnostic output
c capta          : matrix Theta, at "cntr"
c bfctr          : constant factor to compute B from Theta
c capt           : current Theta matrix
c capb           : current B matrix
c vcapba,vcapb   : eigenvalues for capba, capb
c cmplx1         : complex constant 1+0i
c cmplxh         : complex "version" of the stepsize h, h+0i
c h12, h13, h16  : complex  1/2*h, 1/3*h, 1/6*h used for Runge-Kutta
c 
c--------------------------------------------------------------------
c Data initialization
c
c
      cmplx1=dcmplx(1.0d0,0.0d0)
      mesh=1000 
      submsh=8
c
      h=(b-cntr)/dfloat(mesh*(submsh+1))
      cmplxh=dcmplx(h,0.0d0)
      h12=dcmplx(h/2.0d0,0.0d0)
      h16=dcmplx(h/6.0d0,0.0d0)
      h13=dcmplx(h/3.0d0,0.0d0)
c--------------------------------------------------------------------
c
      do 101 i=1,2
        qdold(i)=qda(i)
        do 101 j=1,2
          capt(i,j)=capta(i,j)
  101 continue
      icnt=0
c
c--------------------------------------------------------------------
c --- we move from cntr to b in steps of h
c     at each stage we compute the (approx) new capt 
c     after we have done the submesh points we find capb
c     and find the arguments of the eigenvalues of capb
c     we keep track of the total number of times that they
c     pass across the pos real axis ( in anti-cw dir )
c
      iter=0
      do 500 indx1=1,mesh
        do 430 indx2=0,submsh
          x=cntr+dfloat(iter)*h
          iter=iter+1
          call findg(x,lamda,capt,t1)
          call matadd(cmplx1,capt,h12,t1,w1)
          x1=x+h/2.0d0
          call findg(x1,lamda,w1,t2)
          call matadd(cmplx1,capt,h12,t2,w2)
          call findg(x1,lamda,w2,t3)
          call matadd(cmplx1,capt,cmplxh,t3,w3)
          x2=x+h
          call findg(x2,lamda,w3,t4)
          do 400 i=1,2
            do 400 j=1,2
              capt(i,j)=capt(i,j)+h16*t1(i,j)+h13*t2(i,j)+
     *                   h13*t3(i,j)+h16*t4(i,j)
  400     continue
c
c ---     we are done if this is a submesh point
c
  430   continue
c
c ---   we form capb = bfctr * capt
c
        call matmul(bfctr,0,capt,0,capb)
c
c ---   we get the eigenvalues of capb
c      
        call eigval(capb,vcapb(1),vcapb(2))
c
c ---   we compute the quadrants of the eigenvalues, and compare
c       with the previous values to see if we have crossed the 
c       positive axis
c        
        do 435 i=1,2
          qd(i)=quad(vcapb(i))
  435   continue
c
        do 450 i=1,2
          if((qdold(i).eq.4).and.(qd(i).eq.1)) then  
            icnt=icnt+1
          else
            if((qdold(i).eq.1).and.(qd(i).eq.4)) then   
              icnt=icnt-1
            endif
          endif
  450   continue
c
        do 460 i=1,2
          qdold(i)=qd(i)
  460   continue
c--------------------------------------------------------------------
c --- we are done with the integration, so we send back the results
c
  500 continue
      kount=icnt
      return
      end
c
c********************************************************************
c
c  FINDG
c
c this is a routine to compute the matrix G in the diff eq
c      '
c theta  = G
c
c G= i * theta * omega
c
c********************************************************************
c
      subroutine findg(x,lamda,capt,capg)
c
      double precision x,lamda
      double complex capt(2,2),capg(2,2),
     * tmp1(2,2),tmp2(2,2),tmp3(2,2),
     * tmp4(2,2)
      double complex invp(2,2),capq(2,2),
     * capw(2,2),capr(2,2)
      double complex cmplx1,clam,halfi
c--------------------------------------------------------------------
c Data initialization
c
      cmplx1=dcmplx(1.0d0,0.0d0)
      clam=dcmplx(lamda,0.0d0)
      halfi=dcmplx(0.0d0,0.5d0)
c
c --- we begin by getting the coefficient matrices
c
      call pinv(x,invp)
      call qmat(x,capq)
      call wmat(x,capw)
      call matadd(cmplx1,capq,-clam,capw,capr)
c
c --- we form tmp1=theta + I, tmp2=theta-I
c
      do 110 i=1,2
        do 100 j=1,2
          tmp1(i,j)=capt(i,j)
          tmp2(i,j)=capt(i,j)
  100   continue
        tmp1(i,i)=tmp1(i,i)+cmplx1
        tmp2(i,i)=tmp2(i,i)-cmplx1
  110 continue
c                                   *
c --- form tmp3=invp*tmp1, tmp4=tmp1 *tmp3
c
      call matmul(invp,0,tmp1,0,tmp3)
      call matmul(tmp1,1,tmp3,0,tmp4)
c                                   *
c --- form tmp1=capr*tmp2, tmp3=tmp2 *tmp1
c
      call matmul(capr,0,tmp2,0,tmp1)
      call matmul(tmp2,1,tmp1,0,tmp3)
c
c --- form tmp1= 1/2i*tmp4 - 1/2i*tmp3   this is i*omega
c
      call matadd(halfi,tmp4,-halfi,tmp3,tmp1)
c
c --- form capg=capt*tmp1
c
      call matmul(capt,0,tmp1,0,capg)
c
c --- we are done
c
      return
      end
c
c
c********************************************************************
c
c   QUAD
c
c this is a function that returns the number of the quadrant 
c the complex number is in
c
c********************************************************************
c
      function quad(z)
c
      integer quad
      double precision rez,imz
      double complex z
c
      rez=dble(z)
      imz=dimag(z)
c
      if(rez.ge.0.0d0) then
        if(imz.ge.0.0d0) then
          quad=1
        else
          quad=4
        endif
      else
        if(imz.ge.0.0d0) then
          quad=2
        else
          quad=3
        endif
      endif
      return
      end
c
c
c
c********************************************************************
c
c   MATADD
c
c  this routine forms the scaled sum of two complex 2x2 matrices
c 
c  c = k1*a + k2*b
c
c********************************************************************
c
      subroutine matadd(k1,a,k2,b,c)
c
      double complex k1,k2,a(2,2),b(2,2),c(2,2)
c
      do 100 i=1,2
        do 100 j=1,2
          c(i,j)=k1*a(i,j)+k2*b(i,j)
  100 continue
c
      return
      end
c
c
c********************************************************************
c
c   INVERT
c
c  this routine forms the inverse of a 2x2 complex matrix
c     -1
c  b=a
c
c********************************************************************
c
      subroutine invert(a,b,ifail)
c
      double complex deta,a(2,2),b(2,2),cmplx0
      integer ifail
c
      cmplx0=dcmplx(0.0d0,0.0d0)
c
      deta=a(1,1)*a(2,2)-a(1,2)*a(2,1)
      if (deta.ne.cmplx0) then
        b(1,1)=a(2,2)/deta
        b(1,2)=-a(1,2)/deta
        b(2,1)=-a(2,1)/deta
        b(2,2)=a(1,1)/deta
      else
        ifail=1
      endif
c
      return
      end
c
c
c********************************************************************
c
c  EIGVAL 
c
c  this routine finds the eigenvalues of a 2x2 complex matrix
c
c********************************************************************
c
      subroutine eigval(mata,e1,e2)
c
      double complex mata(2,2),discr,e1,e2,cmplx2,cmplx4,a,b,c,d
c
      cmplx2=dcmplx(2.0d0,0.0d0)
      cmplx4=dcmplx(4.0d0,0.0d0)
      a=mata(1,1)
      b=mata(1,2)
      c=mata(2,1)
      d=mata(2,2)
c
      discr=zsqrt((a+d)*(a+d)-cmplx4*(a*d-b*c))
      e1=(a+d+discr)/cmplx2
      e2=(a+d-discr)/cmplx2
c
      return
      end
c
c
c
c
c********************************************************************
c
c MATMUL
c
c this routine forms the product of two 2x2 complex matrices
c (or their adjoints)
c     if aflg=0 use matrix a, if aflg=1 use the adjoint of a
c     if bflg=0 use matrix b, if bflg=1 use the adjoint of b
c
c********************************************************************
c
      subroutine matmul(a,aflg,b,bflg,c)
c
      integer aflg,bflg
      double complex a(2,2),b(2,2),c(2,2)
c
      if(aflg.eq.0) then
        if(bflg.eq.0) then
          do 100 i=1,2
            do 100 j=1,2
              c(i,j)=dcmplx(0.0d0,0.0d0)
              do 100 k=1,2
                c(i,j)=c(i,j)+a(i,k)*b(k,j)
  100     continue
        else
          do 200 i=1,2
            do 200 j=1,2
              c(i,j)=dcmplx(0.0d0,0.0d0)
              do 200 k=1,2
                c(i,j)=c(i,j)+a(i,k)*dconjg(b(j,k))
  200     continue
        endif
      else
        if(bflg.eq.0) then
          do 300 i=1,2
            do 300 j=1,2
              c(i,j)=dcmplx(0.0d0,0.0d0)
              do 300 k=1,2
                c(i,j)=c(i,j)+dconjg(a(k,i))*b(k,j)
  300     continue
        else
          do 400 i=1,2
            do 400 j=1,2
              c(i,j)=dcmplx(0.0d0,0.0d0)
              do 400 k=1,2
                c(i,j)=c(i,j)+dconjg(a(k,i))*dconjg(b(j,k))
  400     continue
        endif
      endif
      return
      end
c
c
c********************************************************************
c
c SETUP
c
c this subroutine defines the interval and initial conditions
c for the problem, and initializes various constants
c
c********************************************************************
c
      subroutine setup 
c
      double precision pi,a,b,cntr
      double complex b1(2,2),b2(2,2),cmplx0,cmplx1 
      double complex bc1(2,2),bc2(2,2)   
      integer track,flag
c
      common/blk1/pi,a,b,cntr,b1,b2,track
c
c --- set the useful complex constants
c
      cmplx0=dcmplx(0.0d0,0.0d0)
      cmplx1=dcmplx(1.0d0,0.0d0)
c
c --- get the values for the left and right endpoints of the
c     interval that the scalar problem is on:  [a,b]
c     and compute the midpoint of the interval, cntr
c
 100  write(6,*)' '
      write(6,*)'Enter the values for the interval endpoints [a,b]'
      write(6,*)' a='
      read(5,*) a
      write(6,*)' b='
      read(5,*) b
      if(b.le.a) goto 100
      cntr=(a+b)/2.0d0
c--------------------------------------------------------------------
c --- get the values of the constant matrices that define the
c     initial conditions for the scalar problem
c
      write(6,*)' '
      write(6,*)'First boundary condition:'
      write(6,*)'  C11*y(a)+C12*py''(a)+D11*y(b)+D12*py''(b)=0'
      write(6,*)'Please enter complex constants C11, C12, D11, D12'
      write(6,*)'  C11='
      read(5,*) bc1(1,1)
      write(6,*)'  C12='
      read(5,*) bc1(1,2)
      write(6,*)'  D11='
      read(5,*) bc2(1,1)
      write(6,*)'  D12='
      read(5,*) bc2(1,2)
      write(6,*)' '
      write(6,*)'Second boundary condition:'
      write(6,*)'  C21*y(a)+C22*py''(a)+D21*y(b)+D22*py''(b)=0'
      write(6,*)'Please enter complex constants C21, C22, D21, D22'
      write(6,*)'  C21='
      read(5,*) bc1(2,1)
      write(6,*)'  C22='
      read(5,*) bc1(2,2)
      write(6,*)'  D21='
      read(5,*) bc2(2,1)
      write(6,*)'  D22='
      read(5,*) bc2(2,2)
c
c --- check that the user entered boundary conditions for a 
c     valid (selfadjoint) problem
c
      call checkbc(bc1,bc2,flag)
      if(flag.ne.0) then
        write(6,*)' '
        write(6,*)' ERROR IN BOUNDARY CONDITIONS'
        write(6,*)' '
        write(6,*)'Proper B.C. must satisfy four conditions:'
        write(6,*)'          ___   ___           ___   ___'
        write(6,*)' (1)  C11*C12 - C11*C12 = D11*D12 - D11*D12'
        write(6,*)'          ___   ___           ___   ___'
        write(6,*)' (2)  C21*C22 - C21*C22 = D21*D22 - D21*D22'
        write(6,*)'          ___   ___           ___   ___'
        write(6,*)' (3)  C21*C12 - C11*C22 = D21*D12 - D11*D22'
        write(6,*)' '
        write(6,*)' (4)  Matrix (C|D) must be full rank'
        write(6,*)' '
        write(6,*)'The B.C. you have entered fails to satisfy'
        write(6,*)'condition ',flag
        stop
      endif
c
c--------------------------------------------------------------------
c --- Define the matrices for the boundary conditions for the sep.
c     matrix problem
c
      b1(1,1)=bc2(1,1)
      b1(1,2)=bc1(1,1)
      b1(2,1)=bc2(2,1)
      b1(2,2)=bc1(2,1)
c
      b2(1,1)=bc2(1,2)
      b2(1,2)=-bc1(1,2)
      b2(2,1)=bc2(2,2)
      b2(2,2)=-bc1(2,2)
c
      return
      end
c
c
c********************************************************************
c
c   PINV
c
c********************************************************************
c
      subroutine pinv(x,invp)
c
      double complex invp(2,2),cmplx1
      double precision x
c
      double precision pi,a,b,cntr,pvalu
      double complex b1(2,2),b2(2,2) 
      integer track
c
      common/blk1/pi,a,b,cntr,b1,b2,track
c
      cmplx1=dcmplx(1.0d0,0.0d0)
c
c --- form the matrix problem coeff.invP(x) using the coeff. function
c     (user defined) for the scalar problem p(x)
c
      pvalu=coefp(x)
      invp(1,1)=cmplx1/dcmplx(pvalu,0.0d0)
      invp(1,2)=dcmplx(0.0d0,0.0d0)
      invp(2,1)=dcmplx(0.0d0,0.0d0)
      pvalu=coefp(a+b-x)
      invp(2,2)=cmplx1/dcmplx(pvalu,0.0d0)
c
      return
      end
c
c
c********************************************************************
c
c   QMAT
c
c********************************************************************
c
      subroutine qmat(x,capq)
c
      double complex capq(2,2)
      double precision x
c
      double precision pi,a,b,cntr,qvalu
      double complex b1(2,2),b2(2,2) 
      integer track
c
      common/blk1/pi,a,b,cntr,b1,b2,track
c
c --- form the matrix problem coeff. Q(x) using the coeff. function
c     (user defined) for the scalar problem q(x)
c
      qvalu=coefq(x)
      capq(1,1)=dcmplx(qvalu,0.0d0)
      capq(1,2)=dcmplx(0.0d0,0.0d0)
      capq(2,1)=dcmplx(0.0d0,0.0d0)
      qvalu=coefq(a+b-x)
      capq(2,2)=dcmplx(qvalu,0.0d0)
c
      return
      end
c
c
c********************************************************************
c
c   WMAT
c
c********************************************************************
c
      subroutine wmat(x,capw)
c
      double complex capw(2,2)
      double precision x
c
      double precision pi,a,b,cntr,wvalu
      double complex b1(2,2),b2(2,2) 
      integer track
c
      common/blk1/pi,a,b,cntr,b1,b2,track
c
c --- form the matrix problem coeff. W(x) using the coeff. function
c     (user defined) for the scalar problem w(x)
c
      wvalu=coefw(x)
      capw(1,1)=dcmplx(wvalu,0.0d0)
      capw(1,2)=dcmplx(0.0d0,0.0d0)
      capw(2,1)=dcmplx(0.0d0,0.0d0)
      wvalu=coefw(a+b-x)
      capw(2,2)=dcmplx(wvalu,0.0d0)
c
      return
      end
c
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c    USER  DEFINED  SUBROUTINES
c
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c
c
c********************************************************************
c
c   coefp
c
c********************************************************************
c
      function coefp(x)
c
      double precision coefp,x
c
      coefp=2.0d0
c
      return
      end
c
c
c********************************************************************
c
c   coefq
c
c********************************************************************
c
      function coefq(x)
c
      double precision coefq,x
c
      coefq=1.0d0+x*x
c
      return
      end
c
c
c********************************************************************
c
c   coefw
c
c********************************************************************
c
      function coefw(x)
c
      double precision coefw,x
c
      coefw=1.0d0
c
      return
      end
c
\end{verbatim}