program fit
c  Finding delta that gives zero eigenvalue in the potential  V(r) = -1/r+lambda(1-exp(-delta r))/r
      implicit real*8 (a-h,o-z)
      common nq,en,enex,smax,alam,delta,v0,v1
      external func
      print*,'Enter Lambda-guess, dlambda'
      read*,alam0,dd
      eps=1.d-8
      eta=0
      print*,'Enter Number of electrons, Quantum number'
      read*,nel,nq
      ifail=0
      call c05agf(alam0,dd,eps,eta,func,a,b,ifail)
      if(ifail.ne.0)print*,'IFAIL =',ifail
      print10,alam0,(nel-1)/alam0
  10  format(' -----',2f13.8)
  20  format(1x,2f13.8)
      stop      
      end
      
      real*8 function func(alam0)           
      implicit real*8 (a-h,o-z)
      dimension xpoint(4),hmax(2,4),y(9),w(9,9)
      common nq,en,enex,smax,alam,delta,v0,v1
      external coeffn,bdyval,monit,report,fcn,fcn1,gcond
      pi=4*datan(1.d0)
      toldefa=1.d-9
      alam=alam0
c Define Delta as a function of Lambda!
      if(nq.eq.1)delta=1.26422-0.408707*alam+0.0251437*alam**2
c      if(nq.eq.1)delta=1.26461-0.410073*alam+0.0261292*alam**2
      if(nq.eq.2)delta=0.530339-0.215727*alam-0.0557214*alam**2
      if(nq.eq.3)delta=0.486165-0.870471*alam+1.10241*alam**2
     . -0.588113*alam**3      
      tol0=toldefa
      tolg=tol0*10000
      toll=tol0/4
      x0=tol0**(1.d0/4)/3
  60  format(9x,a19,e10.2)
      ind=nq-1
c      en0=-(1-alam)**2/2/nq**2
c      if(abs(en0).lt.1.d-3)en0=-1.d-3
      en0=0
      en=en0
      d=(1-alam)**2/2/nq**2
  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
   2  if(v(x2).gt.en)goto1
      x2=0.9*x2
      if(x2.lt.x0)print*,'Root of V(r) was not found!'
      goto2
   1  if(v(x1).gt.en .or. v(x2).lt.en)print*,'Root was not localised!'
      do 101 iter=1,100
      x=(x1+x2)/2
      vx=v(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)
      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
  12  if(v(x2).gt.en)goto11
      x2=0.9*x2
      if(x2.lt.x0)print*,'2 Root of V(r) was not found!'
      goto12
  11  if(v(x1).gt.en .or. v(x2).lt.en)print*,'2 Root was not localised!'
      do 102 iter=1,100
      x=(x1+x2)/2
      vx=v(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)
      en=elam
      d=delam*20
c      d=delam*2
c   Repeating the calculation with accelerating accuracy
c   Xend is determined from eq. exp(-q.cl. int.)=sqrt(toll)
      smax=-dlog(toll)/2
      x1=x0
      x2=1.d9
  22  if(v(x2).gt.en)goto21
      x2=0.9*x2
      if(x2.lt.x0)print*,'3 Root of V(r) was not found!'
      goto22
  21  if(v(x1).gt.en .or. v(x2).lt.en)print*,'3 Root was not localised!'
      do 103 iter=1,100
      x=(x1+x2)/2
      vx=v(x)
      if(vx.lt.en)x1=x
 103  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=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,'lam,del,xend,en',alam,delta,xend,elam
      en=elam
  40  format(f10.5,f16.8)
  10  format(1x,a15,4f12.6)
      func=en
      return
      end
            subroutine coeffn(p,q,dqdl,x,elam,jint)
      implicit real*8 (a-h,o-z)
      p=-.5d0
      dqdl=-1.d0
      q=v(x)-elam
      return
      end
            subroutine bdyval(xl,xr,elam,yl,yr)
      implicit real*8 (a-h,o-z)
      dimension yl(3),yr(3)
      common nq,en,enex,smax,alam,delta,v0,v1
      h00=1
      h01=v0*h00
      h02=(v0*h01+(v1-en)*h00)/3
      fun=xl*(h00+h01*xl+h02*xl**2)
      dfun=h00+2*h01*xl+3*h02*xl**2
      yl(1)=fun
      yl(2)=dfun/2
      yr(1)=1
      p2=2*(v(xr)-elam)
      dp2=2*vdif(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
c      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 nq,en,enex,smax,alam,delta,v0,v1
      v0=-1
      v1=alam*delta
      v=-1/x+alam*(1-exp(-delta*x))/x
      return
      end
            real*8 function vdif(x)
      implicit real*8 (a-h,o-z)
      common nq,en,enex,smax,alam,delta,v0,v1
      vdif=(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 nq,en,enex,smax,alam,delta,v0,v1
      gcond=y(1)-smax
      return
      end
            subroutine fcn(t,y,f)
      implicit real*8 (a-h,o-z)
      dimension y(3),f(3)
      common nq,en,enex,smax,alam,delta,v0,v1
	f(1)=y(2)
      f(2)=2*(v(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 nq,en,enex,smax,alam,delta,v0,v1
      a=v(t)-en
      f(1)=0
      if(a.le.0)return
      f(1)=dsqrt(2*a)
      return
      end