subroutine boundsch(n,eig,nmx,wfn,l,nz,nroot,rv,h,v0,v0p,
     x qf,emin,eguess,niter,ierr,ntroot,a,b,p1,p2)
c  pgm to solve radial schroedinger equation for nroot bound state
c    energies and wavefunctions for angular momentum l
c    with potential rv/r, given in uniform mesh of n points
c   r=i*h, i=1,...n-1 ;assuming p(r)=C*r**(l+1)*polynomial(r) for r==0; 
c                               p((n+1)*h)=0
c  nz=nuclear charge
c  emin=is estimate of lowest eigenvalue; used if nz=0
c     otherwise, set to the value of -(nz/(l+1))**2
c
c  eguess is an estimate of eigenvalues (used only if consistent
c    with emin and emax)
c
c  It is assumed that the wavefunction has np-l-1 nodes, where
c    np is the principle quantum number-- np=1,2,..nroot
c
c  uses Noumerov algorithm 
c
c  For l=0,1 corrections are needed to approximate wfn(r=0)
c     These depend upon:
c         e0 (current guess of energy eigenvalue)
c         l,nz
c         v(0) == v0 electronic potential at r=0
c         v'(0) == v0p derivative of electronic potential at r=0
c
c  Corrections are also needed for r>n*h, depending on:
c         e0 (current guess of energy eigenvalue
c         the extrapolated value of rv == r * v
c
c ierr=an nroot digit number indicating status of each root
c   a digit of 1 indicates success in converging root
c              2 indicates near success in converging root
c              9 indicates that root not found
c
c first check how many roots expected;
c   =  =  ntroot (returned as argument)
c
      implicit real*8 (a-h,o-z)
      parameter (convre=1.d-11,vlrg=1.d30,inc=10)
      dimension rv(n),eig(nroot),wfn(nmx,nroot),eguess(nroot)
      dimension p1(nmx),p2(nmx)
      dimension a(nmx),b(nmx)
      err=n*(nz*h)**5
      convrez=convre
      if (nz.gt.0) convrez=convre*nz
c     write(6,*) 'expected error = ',err
      ierr=0
      iter=0
      if (n.gt.nmx) then
       write(6,*) 'error in boundsch -- n,nmx = ',n,nmx
       stop
       endif
      angm=l*(l+1)
      do 10 i=1,n
      a(i)=angm/(i*i)+h*rv(i)/i
   10 continue
      do 12 i=1,n
      b(i)=-1.2d0+0.1d0*a(i)
   12 continue
      ssl=0.1d0
      ssd=1.d0
      do 14 i=1,n
   14 a(i)=2.4d0+a(i)
      write(6,*) 'z , l = ',nz,l
c check how many roots expected by integration outward at
c   energy = 0
       energy = 0
c
c  start outward integration
c    correct behavior near r=0
      b0p0=0
      if (l.gt.1) go to 101
      c1=-nz/(l+1.d0)
      c2=((v0-energy)-2*nz*c1)/(4*l+6.d0)
      c3=(v0p+(v0-energy)*c1-2*nz*c2)/(6*l+12.d0)
      fact=1.d0+h*(c1+h*(c2+h*c3))
      if (l.eq.0) then
        b0p0=-0.2d0*h*nz/fact
        endif
      if (l.eq.1) then
        b0p0=0.2d0/fact
        endif
  101 continue
      eh2=0
      node=0
      p1(1)=h**(l+1)
      p1(2)=(a(1)+b0p0-eh2*ssd)*p1(1)/(eh2*ssl-b(2))
       do j=3,n
       p1(j)=((a(j-1)-eh2*ssd)*p1(j-1)+
     x    (b(j-2)-eh2*ssl)*p1(j-2))/(eh2*ssl-b(j))
       if ((j.lt.n-10).and.p1(j)*p1(j-1).le.0.d0) node=node+1
c  if wavefunction too large, renormalize
      if (dabs(p1(j)).gt.vlrg) then
       scale=1.d0/vlrg
       ! call dscal(j,scale,p1,1)
       p1(1:j)=scale*p1(1:j)
       endif
       enddo
       write(6,*) ' nodes at e=0  ', node
       mxroot=node+1
       ntroot=node
       if (mxroot.lt.nroot) then
       write(6,*)'error in boundsch - for l = ',l
       write(6,*) nroot,' states requested but only',mxroot,' possible'
       do ir=mxroot+1,nroot
       ierr=ierr+9*(10**(ir-1))
       enddo
       endif
       mxroot=min0(mxroot,nroot)
c
       if (nz.eq.0) energy=-dabs(emin)
       if (nz.ne.0) energy=-(nz/(l+1.d0))**2
       emin=energy-err
       emax=0.d0
       iroot=1
       best=1.d10
       energy=eguess(iroot)
       if (energy.lt.emin) energy=emin
       if (energy.gt.emax) energy=emax
 3001 continue
c     write(6,*) 'iter,iroot,energy',iter,iroot,energy
c     write(6,*) 'emin,max',emin,emax
      eh2=energy*h*h
c  start inward integration
      ! call dcopy (n,0.d0,0,p2,1)
      p2=0.d0
c  start integration at n
      least=min0(10,n/4)
      r=n*h
      r2=r*r
      veff=rv(n)/n+rv(n-1)/(n-1)
      veff=0.5d0*veff/h+angm/r2
      ppp=dsqrt(dabs(veff-energy))
      j=n
      p2(j)=1.d0/dsqrt(ppp)
       rvp1=2*rv(j)-rv(j-1)
       bnp1=-1.2d0+0.1d0*(h*rvp1/(j+1)+angm/(j+1)**2)
      r=r+h
      r2=r*r
      pppp1=dsqrt(dabs(rvp1/r+angm/r2-energy))
      arg=0.5d0*(ppp+pppp1)*h
      pnp1=dexp(-arg)/dsqrt(pppp1)
      p2(j-1)=((a(j)-eh2*ssd)*p2(j)+
     x   (bnp1-eh2*ssl)*pnp1)/(eh2*ssl-b(j-1))
      j=j-1
 4001 j=j-1
      p2(j)=((a(j+1)-eh2*ssd)*p2(j+1)+
     x   (b(j+2)-eh2*ssl)*p2(j+2))/(eh2*ssl-b(j))
c  if wavefunction too large, renormalize
      if (dabs(p2(j)).gt.vlrg) then
       many=n-j+1
       scale=1.d0/vlrg
       ! call dscal(many,scale,p2(j),1)
        p2(j:n)=scale*p2(j:n)
       endif
c  step inward until p2 stops increasing or j=least
      if (p2(j).gt.p2(j+1).and.j.gt.least) go to 4001
      match=j+1
      rin=(p2(match+1)-p2(match-1))/(2*h*p2(match))
c     write(6,*) ' match point = ',match,rin,p2(match)
c  start outward integration
c    correct behavior near r=0
      b0p0=0
      if (l.gt.1) go to 1001
      c1=-nz/(l+1.d0)
      c2=((v0-energy)-2*nz*c1)/(4*l+6.d0)
      c3=(v0p+(v0-energy)*c1-2*nz*c2)/(6*l+12.d0)
      fact=1.d0+h*(c1+h*(c2+h*c3))
      if (l.eq.0) then
        b0p0=-0.2d0*h*nz/fact
        endif
      if (l.eq.1) then
        b0p0=0.2d0/fact
        endif
 1001 continue
      node=0
      p1(1)=h**(l+1)
      p1(2)=(a(1)+b0p0-eh2*ssd)*p1(1)/(eh2*ssl-b(2))
       do j=3,match+1
       p1(j)=((a(j-1)-eh2*ssd)*p1(j-1)+
     x    (b(j-2)-eh2*ssl)*p1(j-2))/(eh2*ssl-b(j))
       if (p1(j)*p1(j-1).le.0.d0) node=node+1
       enddo
      rout=(p1(match+1)-p1(match-1))/(2*h*p1(match))
c     write(6,*) 'node,match,rin,rout',node,match,rin,rout 
c check whether node = (iroot-1)
c   not enough nodes -- raise energy
      if (node.lt.iroot-1) then
       emin=dmax1(emin,energy)-err
       energy=emax-(emax-energy)*ranx()
       ifac=9
       go to 2001
       endif
c   too many nodes -- lower energy
      if (node.gt.iroot-1) then
       if (energy.le.emin) then
        ierr=ierr+9*(10**(iroot-1))
        write(6,*) 'error in boundsch -- emin too high',l,nz,energy
        endif
       emax=dmin1(emax,energy)+err
       energy=emin+(energy-emin)*ranx()
       go to 3001
       endif
c   correct number of nodes -- estimate correction
       if (node.eq.iroot-1) then
          do j=1,match
          p1(j)=p1(j)/p1(match)
c         write(6,*) 'j,p1',j,p1(j)
          enddo
          do j=match,n
          p1(j)=p2(j)/p2(match)
c         write(6,*) 'j,p2',j,p1(j)
          enddo
          scale=1.d0/overlap(n,h,p1,p1)
          dele=(rout-rin)*scale
c         write(6,*) 'dele,scale',dele,scale
          x=dabs(dele)
          if (x.lt.best) then
          eig(iroot)=energy
          best=x
          endif
        if (dabs(dele).le.convrez) go to 5001
          energy=energy+dele
c if energy is out of range, pick random energy in correct range
          if (emin-energy.gt.convrez.or.energy-emax.gt.convrez)
     x        energy=emin+(emax-emin)*ranx()
           ifac=2
        endif
 2001     iter=iter+1
          if (iter.gt.niter) then
          ierr=ierr+ifac*(10**(iroot-1))
          write(6,*) 'no convergence in boundsch',iroot,l,dele,energy
          write(6,*) ' best guess of eig, dele = ',eig(iroot),best
          if (iroot.lt.mxroot) then
          do ir=iroot+1,mxroot
          ierr=ierr+9*(10**(ir-1))
          enddo
          endif
c         write(6,*) 'returning from boundsch -- ierr=',ierr
          return
          endif
          go to 3001
 5001 continue
c  eigenvalue found
      eig(iroot)=energy
      scale=dsqrt(scale)
      ! call dscal(n,scale,p1,1)
       p1=scale*p1
      ! call dcopy(n,p1,1,wfn(1,iroot),1)
       wfn(:,iroot)=p1
c     write(6,*) 'root found',iroot,energy,scale
      ierr=ierr+10**(iroot-1)
      iroot=iroot+1
      if (iroot.le.mxroot) then
      emin=energy-err
      emax=0
      energy=eguess(iroot)
       if (energy.lt.emin) energy=emin
       if (energy.gt.emax) energy=emax
      best=1.d10
      go to 3001
      endif
c  all roots found -- return to calling program
c         write(6,*) 'returning from boundsch -- ierr=',ierr
      return
      end