program screened
c   Arbitrary angular momentum      
      implicit real*8 (a-h,o-z)
      dimension xpoint(4),hmax(2,4),y(9),w(9,9)
      common en,smax,alam,delta,v0,v1,lq,qll
      external coeffn,bdyval,monit,report,fcn,fcn1,gcond
  50  format(' Solving the Schrodinger equation for bound states'/
     .'in a potential V(r) = -1/r+lambda(1-exp(-delta r))/r')
      pi=4*datan(1.d0)
      toldefa=1.d-9
      print*,'Files: 1 - output1 (acc.=tol), 2 - output2 (acc.=tol/4)'
c      open(1,file='temp1.out')
c      open(2,file='temp2.out')
c      open(3,file='screen.dat')
      write(1,50)
      write(2,50)
      print*,'Type: n, l'
      read*,nq,lq
      qll=lq*(lq+1)/2.d0
      print*,'Type: Delta'
      read*,delta
      print*,'Type accuracy, if 0 then',toldefa
      tol0=toldefa
      read*,tol0
      if(tol0.eq.0)tol0=toldefa
      tolg=tol0*10000
      toll=tol0/4
      write(1,60)'Accuracy (tol) =',tol0
      write(2,60)'Accuracy (tol/4) =',toll
      x0=tol0**(1.d0/4)/3
  60  format(9x,a19,e10.2)
      ind=nq-lq-1
      print*,'Type: Lambda - min, max, step'
      read*,alamin,alamax,alastep
      en=0
      do 1 alam=alamin,alamax,alastep
c !      
      x1=1.d0/2
      x2=1.d0/3
      y1=0.596071
      y2=0.614559
      delta=(y1*(alam-x2)-y2*(alam-x1))/(x1-x2)
c      en0=-(1-alam)**2/2/nq**2
c !
      en0=en
       if(abs(en0).lt.1.d-3)en0=-1.d-3
      en=en0
c      d=(1-alam)**2/2/nq**2
c !
      d=1.d-3
      print30,alam,nq,en,d
      print30,alam,nq,en,d
  30  format(' lambda, q.number, Ecl, dE',f6.3,i4,2f20.12)
c   Xend is determined from eq. exp(-q.cl. int.)=sqrt(tolg)
      smax=-dlog(tolg)/2
      x1=x0
      x2=1.d9
      if(veff(x1).gt.en .or. veff(x2).lt.en)
     .  print*,'Root was not localised!'
      do 101 iter=1,100
      x=(x1+x2)/2
      vx=veff(x)
      if(vx.lt.en)x1=x
 101  if(vx.ge.en)x2=x
      y(1)=0
      tol=tolg
      ifail=0
      xend=10000
      call d02bhf(x,xend,1,y,tol,0,0.d0,fcn1,gcond,w,ifail)
      xend=x
      xpoint(1)=x0
      xpoint(2)=x0
      xpoint(3)=xend
      xpoint(4)=xend
      match=0
      tol=tolg
      elam=en-d
      delam=2*d
      hmax(1,1)=0
      maxit=100
      maxfun=100 000
      ifail=0
      call d02kef(xpoint,4,match,coeffn,bdyval,ind,tol,elam,delam,
     .        hmax,maxit,maxfun,monit,report,ifail)
      print10,'xend, en0, den =',xend,elam,delam
      en=elam
      d=delam*2
c   Xend is determined from eq. exp(-q.cl. int.)=sqrt(tol0)
      smax=-dlog(tol0)/2
      x1=x0
      x2=1.d9
      if(veff(x1).gt.en .or. veff(x2).lt.en)
     .   print*,'Root was not localised!'
      do 102 iter=1,100
      x=(x1+x2)/2
      vx=veff(x)
      if(vx.lt.en)x1=x
 102  if(vx.ge.en)x2=x
      y(1)=0
      tol=tolg
      ifail=0
      xend=10000
      call d02bhf(x,xend,1,y,tol,0,0.d0,fcn1,gcond,w,ifail)
      xend=x
      xpoint(1)=x0
      xpoint(2)=x0
      xpoint(3)=xend
      xpoint(4)=xend
      match=0
      tol=tol0
      elam=en-d
      delam=2*d
      hmax(1,1)=0
      maxit=100
      maxfun=100 000
      ifail=0
      call d02kef(xpoint,4,match,coeffn,bdyval,ind,tol,elam,delam,
     .        hmax,maxit,maxfun,monit,report,ifail)
      print10,'xend, en, den =',xend,elam,delam
      en=elam
      d=delam*20
c      d=delam*2
      write(1,20)alam,nq,en0,en
  20  format(' lambda=',f6.3,'  n=',i3,'   Ecl= [',f9.4,']   E=',
     . f17.12,8h   psi'=,f16.12)
c   Repeating the calculation with accelerating accuracy
c   Xend is determined from eq. exp(-q.cl. int.)=sqrt(toll)
c      smax=-dlog(toll)/2
c      x1=x0
c      x2=1.d9
c      if(veff(x1).gt.en .or. veff(x2).lt.en)
c     .   print*,'Root was not localised!'
c      do 103 iter=1,100
c      x=(x1+x2)/2
c      vx=v(x)
c      if(vx.lt.en)x1=x
c 103  if(vx.ge.en)x2=x
c      y(1)=0
c      tol=tolg
c      ifail=0
c      xend=10000
c      call d02bhf(x,xend,1,y,tol,0,0.d0,fcn1,gcond,w,ifail)
c      xend=x
      xpoint(1)=x0
      xpoint(2)=x0
      xpoint(3)=xend
      xpoint(4)=xend
      match=0
      tol=toll
      elam=en-d
      delam=2*d
      hmax(1,1)=0
      maxit=100
      maxfun=100 000
      ifail=0
      call d02kef(xpoint,4,match,coeffn,bdyval,ind,tol,elam,delam,
     .        hmax,maxit,maxfun,monit,report,ifail)
      print10,'xend, en, den =',xend,elam,delam
      en=elam
      write(2,20)alam,nq,en0,en
      write(3,40)alam,-en
c      write(3,40)alam,(en+(1-alam)**2/2)/(1-alam)**3
  40  format(f10.5,f16.8)
  10  format(1x,a15,4f19.13)
   1  continue
      end
            subroutine coeffn(p,q,dqdl,x,elam,jint)
      implicit real*8 (a-h,o-z)
      common en,smax,alam,delta,v0,v1,lq,qll
      p=-.5d0
      dqdl=-1.d0
      q=veff(x)-elam
      return
      end
            subroutine bdyval(xl,xr,elam,yl,yr)
      implicit real*8 (a-h,o-z)
      dimension yl(3),yr(3)
      common en,smax,alam,delta,v0,v1,lq,qll
      h00=1
      h01=v0*h00/(lq+1)
      h02=(v0*h01+(v1-en)*h00)/(2*lq+3)
      fun=xl**(lq+1)*(h00+h01*xl+h02*xl**2)
      dfun=xl**lq*((lq+1)*h00
     .    +(lq+2)*h01*xl+(lq+3)*h02*xl**2)
      yl(1)=fun
      yl(2)=dfun/2
      yr(1)=1
      p2=2*(veff(xr)-elam)
      dp2=2*veffdif(xr)
      y=p2
      if(y.lt.0)y=-y
      p=dsqrt(y)
      yr(2)=-(dp2/p2/4+p)
      return
      end
            subroutine monit(maxit,iflag,elam,finfo)
      implicit real*8 (a-h,o-z)
      dimension finfo(15)
      nfun=finfo(4)+.1
      print10,maxit,nfun,elam,iflag
  10  format('+',29x,'it,nfun,elam,ifl=',i3,i7,f20.12,i3)
      return
      end
            subroutine report(x,vv,jint)
      implicit real*8 (a-h,o-z)
      dimension vv(3)
      return
      end
            real*8 function v(x)
      implicit real*8 (a-h,o-z)
      common en,smax,alam,delta,v0,v1,lq,qll
      v0=-1
      v1=alam*delta
      v=-1/x+alam*(1-exp(-delta*x))/x
      return
      end
            real*8 function veff(x)
      implicit real*8 (a-h,o-z)
      common en,smax,alam,delta,v0,v1,lq,qll
      v0=-1
      v1=alam*delta
      veff=qll/x**2-1/x+alam*(1-exp(-delta*x))/x
      return
      end
            real*8 function vdif(x)
      implicit real*8 (a-h,o-z)
      common en,smax,alam,delta,v0,v1,lq,qll
      vdif=(1-alam)/x**2+alam*exp(-delta*x)/x**2*(1+delta*x)
      return
      end
            real*8 function veffdif(x)
      implicit real*8 (a-h,o-z)
      common en,smax,alam,delta,v0,v1,lq,qll
      veffdif=-2*qll/x**3+(1-alam)/x**2
     .   +alam*exp(-delta*x)/x**2*(1+delta*x)
      return
      end
            real*8 function gcond(t,y)
      implicit real*8 (a-h,o-z)
      dimension y(1)
      common en,smax,alam,delta,v0,v1,lq,qll
      gcond=y(1)-smax
      return
      end
            subroutine fcn(t,y,f)
      implicit real*8 (a-h,o-z)
      dimension y(3),f(3)
      common en,smax,alam,delta,v0,v1,lq,qll
	f(1)=y(2)
      f(2)=2*(veff(t)-en)*y(1)
      f(3)=y(1)**2
      return
      end
            subroutine fcn1(t,y,f)
      implicit real*8 (a-h,o-z)
      dimension y(1),f(1)
      common en,smax,alam,delta,v0,v1,lq,qll
      a=veff(t)-en
      f(1)=0
      if(a.le.0)return
      f(1)=dsqrt(2*a)
      return
      end