c *********************************************************************
      subroutine scheq(zz,en,lambda,nofl,kkk,mess,scf,thresh,iflag,npr
     1,de)
      implicit double precision (a-h,o-z)
c     subroutine scheq
c     compute energy eigenvalue and wave function
c      originally written by sherwood skillman
c      rca laboratories, princeton, new jersey, spring 1961
c      modified by frank herman, summer 1961
c      further modified by richard kortum and paul kelly, lockheed
c       research laboratories, palo alto, california,  summer 1962
c further modified by r.m. nieminen, otaniemi,spring 1978
c and last but not least modified by m.p. winter 1979
      dimension qq(561),p(5),q(5),t(5),d(5)
      common/sc/gr(561),r(561),snlo(561),nbl
      common/pot/v(561)
      imp=1
      nb1=nbl+1
      nmb=nb1/5
      if (de.lt.1.d-10)imp=0
      if(imp.eq.1)de=dmax1(de,-e/500.d0)
      eg=0.d0
       z=zz
       lam=lambda
       nn=nofl
       mesh=mess
       c=scf
          i0=0
      many = 150
      iflag=0
   73 e=en
      morev=0
      lessv=0
      emore=0.0d0
      eless=0.0d0
      more=0
      less=0
      mp=0
      nprint=0
      lamm=lam-1
      lamp=lam+1
      xlp=lamp
      ndcr=nn-lamp
      b=lam*lamp
       oc=r(2)
      h=oc
      hsq=h*h
      b3=(v(3)-v(2))/h-z/hsq
      y=h+h
      flps=4*lam+6
      slpt=6*lam+12
      elpt=8*lam+20
      a1=-z/xlp
      ysq=y*y
       b1=-z-z
      ab1=a1*b1
      ab3=a1*b3
c     raise h and y to lam+1
      htl=h
      ytl=y
      if(lam)77,1102,1100
 1100 do 1101 i=1,lam
      htl=htl*h
 1101  ytl=ytl*y
 1102 h1=hsq
      bohs=b/hsq
      boh=b1/h
      bth=b3*h
      bq3=bohs+boh+bth
      bq4=bohs/4.d0+boh/2.d0+bth+bth
      epl=8+lam
      fpl=5+lam
      xifc=c*.21701389d-4
c     start outward integration
   10 nprint=nprint+1
      eps =e-eg
      eg =e
      if(many-nprint)  900,75,75
900   if(mp.eq.2)go to 300
      write (6,901)nn,lam ,z
      call flush(6)
  901 format (21h    no convergence on,i4,i1,f4.0)
      iflag=1
      en=0.d0
      npr=nprint
      return
  77  nstop=7
 777  write(66,770) nstop
 770  format(' error in scheq',i6)
      en=0.d0
      iflag=1
      npr=nprint
      return
   75 do 11 i=1,mesh
   11 snlo(i)=0.0d0
      if(nprint-1)  77,500,529
  500 continue
      do 502  i= 4,mesh
      qq(i) = v(i)+b/(r(i)*r(i))-e
  502 continue
  516 m= mesh
      do 520 i=4,mesh
      if(qq(m)) 519,520,520
  519 ik=m+1
      go to 525
  520 m=m-1
      write(66,562)lam,qq(m+1),e
562   format(/'  l: ',i5,'  q:  ',g15.6,'  e:  ',g15.6/)
  521 nstop =521
c     q is everywhere positive
      go to 777
  525 if(mesh-ik) 526,526,535
  526 eps = qq(mesh-nbl)
      e = e+eps
  529 continue
      do 530 i=4,mesh
  530 qq(i) = qq(i)-eps
      go to 516
  535 continue
   14 ncross=0
      sig=1.0d0
      h=oc
      y=h+h
c     b= lam*(lam+1)
c     b1= -2.d0*z
      b2=3.0d0*z/h-e+2.0d0*v(2)-v(3)
c     b3=(v(3)-v(2))/h -z/hsq
c     a1= -z/(lam+1)
      a2=(ab1+b2)/flps
c     a2=(a1*b1+b2)/(4*lam+6)
      a3=(a2*b1+a1*b2+b3)/slpt
c     a3=(a2*b1+a1*b2+b3)/(6*lam+12)
      a4=(a3*b1+a2*b2+ab3)/elpt
c     a4=(a3*b1+a2*b2+a1*b3)/(8*lam+20)
      p(3)=(1.0d0+h*(a1+h*(a2+h*(a3+h*a4))))*htl
c     p(3)=(1.0d0+a1*h+a2*h**2+a3*h**3+a4*h**4)*h**(xl+1.0d0)
      p(4)=(1.0d0+y*(a1+y*(a2+y*(a3+y*a4))))*ytl
c     p(4)=(1.0d0+a1*y+a2*y**2+a3*y**3+a4*y**4)*y**(xl+1.0d0)
      q(3)=bq3+b2
c     q(3)=(b+b1*h+b2*h**2+b3*h**3)/h**2
      q(4)=bq4+b2
c     q(4)=(b+b1*y+b2*y**2+b3*y**3)/y**2
      snlo(2)=p(3)
      snlo(3)=p(4)
      i=3
       dx=oc
      h1=h**2
      h2=h1/12.0d0
      t(3)=p(3)*(1.0d0-h2*q(3))
      t(4)=p(4)*(1.0d0-h2*q(4))
      d(4)=t(4)-t(3)
      ncount=3
      nint=2
   15 i=i+1
c     if end of mesh is reached, modify trial eigenvalue
      if(i-mesh) 16,3,3
    3 if(ndcr-ncross) 32,33,33
c     return to begininning of outward integration if necessary
   16 q(5) =qq(i)
      if(ik-i)  29,29,21
   21 d(5)=d(4)+h1*q(4)*p(4)
      t(5)=d(5)+t(4)
      if(1.0d0-abs (h2*q(5))) 3,3,501
  501 p(5)=t(5)/(1.0d0-h2*q(5))
      snlo(i) = p(5)
       if(sig) 211,77,212
  211  if(p(5))  23,23,22
  212  if(p(5)) 22,23,23
   22 ncross=ncross+1
c     count changes in sign
      sig=-sig
   23 ncount=ncount+1
      if(7-ncount)77,24,25
   24 ncount=2
   25 nint=nint+1
      if(nbl-nint)77,26,27
   26 dx=dx+dx
      h=dx
      h1=h**2
      h2=h1/12.0d0
      nint=0
      t(5)=p(5)*(1.0d0-h2*q(5))
      t(3)=p(3)*(1.0d0-h2*q(3))
      d(5)=t(5)-t(3)
   27 do 28 k=1,4
      p(k)=p(k+1)
      t(k)=t(k+1)
      d(k)=d(k+1)
   28 q(k)=q(k+1)
      go to 15
   29 if(ncount-2)77,30,21
   30 if(nint-4)21,21,31
c     matching radius has been reached going out
c     if ndcr not equal to ncross,modify trial eigenvalue
   31 eigen=e
      if(ndcr-ncross) 32,35,33
   32 more=1
      mp=0
c     too many crossings, increase absf(e)
      morev=morev+1
       if(morev-1) 50,53,52
   50 nstop=50
      go to 777
   52 if(e -emore) 53,54,54
   53 emore=e
   54  if (less) 55,56,64
   55 nstop=55
      go to 777
   56 e=1.25d0*eg
      go to 10
   33 less=1
      mp=0
c     too few  crossings, decrease abs (e)
      lessv=lessv+1
       if(lessv-1) 57,60,59
   57 nstop=57
      go to 777
   59 if(eless- e) 60,61,61
   60 eless=e
   61 if(more) 62,63,64
   62 nstop=62
      go to 777
   63 e=0.75d0*eg
      go to 10
   64 e=0.5d0*(emore+eless)
      go to 10
   35  if(abs (snlo(i-1))-abs (snlo(i-2))) 351,354,354
c     check to see that wave is in the damped region (absolute value
c     decreasing and signs alike)
  351 if (p(5)) 352,21 ,353
  352 if(snlo(i-2)) 401,21,21
  353 if(snlo(i-2)) 21,21,401
  354 if(1.0d+25 -abs (p(5)))33,33,21
c     large absolute value of p in what should be the damped region
c       indicates too few peaks, decrease absf(e)
c     now ndcr = ncross and matching radius lies in damped region
  401 imatch=i-2
      xmatch=r(i-2)
      ppout=(t(4)-t(2)-0.5d0*(p(4)-p(2)))/h
      s2=ppout/p(3)
c     integration is by 8 applications of newton-cotes closed
c       quadrature for five intervals on each block
c       xifc =(5*h(block 1)/288)/2 ,h(1) =0.0025*scale factor
      sum1=0.0d0
      xif=xifc
      i=1
      value=0.0d0
   36 mm=nmb
      sum2=0.0d0
      xif=xif+xif
   37 y=value
      value=snlo(i+5)**2
      sum2=sum2+19.d0*(value+y)+75.d0*(snlo(i+4)**2+snlo(i+1)**2)
     1 +50.0d0*(snlo(i+2)**2+snlo(i+3)**2)
      i=i+5
       if (imatch-i)  77,39,371
  371  mm=mm-1
      if(mm)77,38,37
   38 sum1=sum2*xif+sum1
       go to 36
   39 sum1= sum1+sum2*xif
   40 s1=sum1/p(3)**2
      pmatch=p(3)
      if(nn-1)77,41,42
   41 xinw=epl* xmatch
c     for n =1, start inward integration at(8+lam)*xmatch or x max
      go to 421
   42 xinw=fpl*xmatch
c     for n not=1, start at (5+lam)*xmatch or x max (end of mesh)
  421 do 44  i= nb1,mesh,nbl
      if(xinw-r(i))  43,43,44
   43 kkk =i
      go to 45
   44 continue
      kkk =mesh
   45 i =kkk
      dx =r(i-1)-r(i)
      h =dx
      xif=0.17361111d-1*dx
      hsq=h*h
      hsq12=hsq/12.0d0
      q(3)= qq(i)
      p(3)= exp (-r(i)*sqrt (q(3)))
  402 sum3=p(3)/q(3)
      i=i-1
      q(4) =qq(i)
  404 p(4)=exp (-r(i)*sqrt (q(4)))
      if (abs (p(4))-1.0d-35) 4041,4041,405
 4041 kkk=kkk -nbl
      if(kkk -imatch) 4042,4042,45
 4042 write (5,4043)z ,nn,lam,kkk
 4043 format  (6hoat z=,f6.0,6h  nl =,i3,i1,7h  kkk =,i5,22h is less tha
     1n imatch =,i5,44h inward integration will be tried at kkk+nbl)
      kkk =kkk +nbl
      p(4) = 1.5d-35
      p(3) = 1.0d-35
  405 if(pmatch)102,77,103
  102 p(3)=-p(3)
      p(4)=-p(4)
  103 snlo(i+1)=p(3)
      snlo(i)=p(4)
      t(3)=p(3)*(1.0d0-hsq12 *q(3))
      t(4)=p(4)*(1.0d0-hsq12*q(4))
      d(4)=t(4)-t(3)
  104 do 106 m=2,nbl
      i=i-1
      q(5) =qq(i)
      d(5)=hsq*q(4)*p(4)+d(4)
      t(5)=d(5)+t(4)
      p(5)=t(5)/(1.0d0-hsq12*q(5))
      if(i-imatch+1)77,200,105
  105 snlo(i)=p(5)
      do 106 k=1,4
      p(k)=p(k+1)
      t(k)=t(k+1)
      d(k)=d(k+1)
  106 q(k)=q(k+1)
      q(5) =qq(i-2)
      d(5)=hsq*q(4)*p(4)+d(4)
      t(5)=d(5)+t(4)
      p(5)=t(5)/(1.0d0-hsq12*q(5))
      p(5)=1.09375d0*p(4)+0.2734375d0*p(5)-0.546875d0*p(3)+0.21875d0*
     1 p(2)-0.0390625d0*p(1)
      i=i-1
      dx=dx/2.0d0
      q(5) =qq(i)
      h=dx
      hsq=h*h
      hsq12=hsq/12.0d0
      t(5)=p(5)*(1.0d0-hsq12*q(5))
      t(4)=p(4)*(1.0d0-hsq12*q(4))
      d(5)=t(5)-t(4)
      snlo(i)=p(5)
      do 107 l=1,4
      p(l)=p(l+1)
      t(l)=t(l+1)
      d(l)=d(l+1)
  107 q(l)=q(l+1)
      go to 104
c     matching radius has been reached coming in
  200 k=kkk
      value=snlo(k)**2
      go to 202
 2001  continue
  201 sum3=sum3+xif*sum4
      xif =xif*0.5d0
  202 mm=nmb
      sum4 =0.0d0
  203 y=value
      value=snlo(k-5)**2
      sum4=sum4+19.d0*(value+y)+75.d0*(snlo(k-1)**2+snlo(k-4)**2)
     1+50.0d0*(snlo(k-2)**2+snlo(k-3)**2)
      k=k-5
       if(k-imatch) 77,2031,2030
 2030  mm=mm-1
       if(mm)  77,2001,203
 2031 sum3=sum3+xif*sum4
  204 s3=sum3/p(4)**2
      ppin=(t(5)-t(3)-0.5d0*(p(5)-p(3)))/h
      s4=ppin/p(4)
      fe=s4-s2
      if(abs(fe).lt.1.d-08)go to 300
c         write(66,4711)mp,nprint,e,fe
          ml=more*less
4711  format(' mp',i2,' iter: ',i4,' e: ',e12.6,' fe: ',e12.6)
      if(mp.eq.2)go to 9120
      if(mp.eq.1)go to 9100
      mp=1
      if(imp.eq.0)de=-e/10.d0
          de=dmin1(de,-e/10.d0)
          if(i0.gt.0)de=de/5**i0
          i0=i0+1
        if(ml.gt.0)de=abs(emore-eless)/10.d0
      if(e.gt.-1.d-35)de=1.d-7
      eold=e
      feold=fe
      if(fe.gt.0)e=e-de
      if(fe.lt.0)e=e+de
      go to 10
9100  if(fe*feold.lt.0)go to 9110
      if(fe.gt.0.d0.and.imp.eq.0) de=-e/10.d0
          if(i0.gt.1.and.ml.gt.0)de=abs(emore-eless)/10.d0
      eold=e
      feold=fe
9101  if(fe.lt.0d0)e=e+de
      if(fe.gt.0d0)e=e-de
      go to 10
9110  continue
9111  mp=2
      fe1=feold
      fe2=fe
      e1=eold
      e2=e
      eold=e
      go to 9130
9120  if(fe1*fe.lt.0d0)e2=e
      if(fe1*fe.lt.0d0)fe2=fe
      if(fe2*fe.lt.0d0)e1=e
      if(fe2*fe.lt.0d0)fe1=fe
      eold=e
6174  format(' e1: ',e12.6,' fe1: ',e12.6,' e2: ',e12.6,
     j' fe2: ',e12.6)
9130  e=(e1*fe2-e2*fe1)/(fe2-fe1)
      if(abs((e2-e1)/e).lt.thresh)go to 300
      if(abs((e-eold)/e).lt.1d-06)go  to 300
      go to 10
  300  pop=pmatch/p(4)
      do 302 j=imatch,kkk
  302 snlo(j)=snlo(j)*pop
      sum1=0.0d0
      j=1
      xif=xifc
      value=0.0d0
  303 mm=nmb
      xif=xif+xif
      sum2=0.0d0
  304 y=value
      value=snlo(j+5)**2
      sum2=sum2+19.d0*(value+y)+75.d0*(snlo(j+4)**2+snlo(j+1)**2)
     1+50.d0*(snlo(j+2)**2+snlo(j+3)**2)
      j=j+5
      mm=mm-1
      if(mm)77,305,304
  305 sum1=sum1+xif*sum2
      if(kkk-j)77,307,303
  307 c1=sqrt (sum1)
      if(snlo(3))308,77,310
  308 c1=-c1
c
c include tail
c
 310  tail=snlo(mesh)**2/(2.d0*sqrt(abs(e)))
      c1=c1+tail
      do 311 i=1,kkk
  311 snlo(i)=snlo(i)/c1
      en=e
      if(abs(en).lt.1.d-6) iflag=1
      npr=nprint
      return
      end
c **********************************************************************
      subroutine simpsh(ff,s)
      implicit double precision (a-h,o-z)
c integrates the matrix ff on the herman-skillman mesh
      common/mess/wj(301),dx,nblock,jblock
      common/sc/gr(561),r(561),snlo(561),nbl
      dimension ff(561)
      nb1=nbl+1
      s=0.d0
      deltax=dx
      do 10 j=1,jblock
      do 11 jk=1,nb1
      i=(j-1)*nbl+jk
  11  s=s+ff(i)*wj(jk)*deltax
  10  deltax=2.d0*deltax
      return
      end
c **********************************************************************
      subroutine scrhs(roo,u,z,q,n)
      implicit double precision (a-h,o-z)
c evaluates the screened poisson integral
c in the herman-skillman mesh
      dimension roo(561),u(561),ss(561),apu(561)
      common/sc/gr(561),r(561),snlo(561),nbl
      pi=4.d0*atan(1.d0)
      c=4.d0*pi
      do 11 i=1,n
  11  ss(i)=-c*roo(i)*r(i)-q**2*u(i)
      do 20 i=1,n
      do 30 j=1,n
      aa=gr(i)/gr(j)
      if(j.gt.i) aa=1.d0/aa
      ab=gr(i)*gr(j)
 30   apu(j)=(aa-ab)*ss(j)
      call simpsh(apu,s)
 20   u(i)=-s/(2.d0*q)-z*gr(i)
      return
      end
c **********************************************************************
      function exc(r,s,ixc)
      implicit double precision (a-h,o-z)
c          ceperley alder perdew zunger
c          or von barth hedin
      data gp,bp1,bp2,gf,bf1,bf2/-.1423d0,1.0529d0,.3334d0,-.0843d0,
     j1.3981d0,.2611d0/
      data ap,bp,cp,dp,af,bf,cf,df/.0311d0,-.048d0,.002d0,-.0116d0,
     j.01555d0,-.0269d0,.0007d0,-.0048d0/
          if(r.lt.1.d-25)go to 100
      pi=4.d0*atan(1.d0)
       if(s.gt.0.99999999d0)s=0.99999999d0
       if(s.lt.-0.99999999d0)s=-0.99999999d0
          sinv=(4.d0*pi*r/3.d0)**(1.d0/3.d0)
          rs=1.d0/sinv
          if(ixc.eq.1)go to 200
      if(rs.lt.1.d0)go to 10
      srs=sqrt(rs)
      excp=gp/(1.d0+bp1*srs+bp2*rs)
      excf=gf/(1.d0+bf1*srs+bf2*rs)
      go to 20
10    aa=log(rs)
      excp=ap*aa+bp+cp*rs*aa+dp*rs
      excf=af*aa+bf+cf*rs*aa+df*rs
20    f=((1.d0+s)**(4.d0/3.d0)+(1.d0-s)**(4.d0/3.d0)-2.d0)/
     1  (2.d0**(4.d0/3.d0)-2.d0)
      excf=excf-.5772521d0/rs
      excp=excp-.4581653d0/rs
      exc=excp+f*(excf-excp)
      return
100   exc=0.d0
      return
200   x=s/2.d0+0.5d0
      d43=4.d0/3.d0
      rsf=rs/75.d0
      rsf2=rsf*rsf
      rsf3=rsf2*rsf
      rsp=rs/30.d0
      rsp2=rsp*rsp
      rsp3=rsp2*rsp
      fcf=(1.d0+rsf3)*log(1.d0+1.d0/rsf)+0.5d0*rsf-rsf2-1.d0/3.d0
      fcp=(1.d0+rsp3)*log(1.d0+1.d0/rsp)+0.5d0*rsp-rsp2-1.d0/3.d0
      epscp=-.0504d0*fcp
      epscf=-.0254d0*fcf
      epsxp=-.91633059d0/rs
      cny=5.1297628d0*(epscf-epscp)
      aa=.5d0**(1.d0/3.d0)
      if(x.lt..000001d0) x=.000001d0
      if(x.gt..999999d0) x=.999999d0
      fx=(x**d43+(1.d0-x)**d43-aa)/(1.d0-aa)
      exc=epsxp+epscp+fx*(cny+4.d0/3.d0*epsxp)/5.1297628d0
      exc=exc/2.d0
      return
      end
c **********************************************************************
      function uxc(r,s,ispin,ixc)
      implicit double precision (a-h,o-z)
c          ceperley alder perdew zunger
      data gp,bp1,bp2,gf,bf1,bf2/-.1423d0,1.0529d0,.3334d0,-.0843d0,
     j1.3981d0,.2611d0/
      data ap,bp,cp,dp,af,bf,cf,df/.0311d0,-.048d0,.002d0,-.0116d0,
     j.01555d0,-.0269d0,.0007d0,-.0048d0/
          if(r.lt.1.d-35)go to 100
      pi=4.d0*atan(1.d0)
          sinv=(4.d0*pi*r/3.d0)**(1.d0/3.d0)
          rs=1.d0/sinv
          if(ixc.eq.1)go to 200
      if(rs.lt.1.d0)go to 10
      srs=sqrt(rs)
      excp=gp/(1.d0+bp1*srs+bp2*rs)
      excf=gf/(1.d0+bf1*srs+bf2*rs)
      uxcp=excp*(1.d0+7.d0/6.d0*bp1*srs+4.d0/3.d0*bp2*rs)
     j/(1.d0+bp1*srs+bp2*rs)
      uxcf=excf*(1.d0+7.d0/6.d0*bf1*srs+4.d0/3.d0*bf2*rs)
     j/(1.d0+bf1*srs+bf2*rs)
      go to 20
10    aa=log(rs)
      excp=ap*aa+bp+cp*rs*aa+dp*rs
      excf=af*aa+bf+cf*rs*aa+df*rs
      uxcp=ap*aa+(bp-ap/3.d0)+2.d0/3.d0*cp*rs*aa+(2.d0*dp-cp)*rs/3.d0
      uxcf=af*aa+(bf-af/3.d0)+2.d0/3.d0*cf*rs*aa+(2.d0*df-cf)*rs/3.d0
20    f=((1.d0+s)**(4.d0/3.d0)+(1.d0-s)**(4.d0/3.d0)-2.d0)/
     j (2.d0**(4.d0/3.d0)-2.d0)
      ddf=4.d0/3.d0*((1.d0+s)**(1.d0/3.d0)-(1.d0-s)**(1.d0/3.d0))
     j /(2.d0**(4.d0/3.d0)-2.d0)
      uxc=uxcp+f*(uxcf-uxcp)+(excf-excp)*(3.d0-2.d0*ispin-s)*ddf
     j -.6108871d0/rs*(1.d0+(3.d0-2.d0*ispin)*s)**(1.d0/3.d0)
      return
100   uxc=0.d0
      return
200   x=s/2.d0+0.5d0
      d43=4.d0/3.d0
      xx=0.5d0-s/2.d0
      rsf=rs/75.d0
      rsf2=rsf*rsf
      rsf3=rsf2*rsf
      rsp=rs/30.d0
      rsp2=rsp*rsp
      rsp3=rsp2*rsp
      fcf=(1.d0+rsf3)*log(1.d0+1.d0/rsf)+0.5d0*rsf-rsf2-1.d0/3.d0
      fcp=(1.d0+rsp3)*log(1.d0+1.d0/rsp)+0.5d0*rsp-rsp2-1.d0/3.d0
      epscp=-.0504d0*fcp
      epscf=-.0254d0*fcf
      epsxp=-.91633059d0/rs
      cny=5.1297628d0*(epscf-epscp)
      aa=.5d0**(1.d0/3.d0)
      if(x.lt..000001d0) x=.000001d0
      if(xx.lt..000001d0) xx=.000001d0
      if(x.gt..999999d0) x=.999999d0
      if(xx.gt..999999d0) xx=.999999d0
      ars=-1.22177412d0/rs+cny
      brs=-0.0504d0*log(1.d0+30.d0/rs)-cny
      trx1=(2.d0*x)**(1.d0/3.d0)
      trx2=(2.d0*xx)**(1.d0/3.d0)
      if(ispin.eq.1)vxc=ars*trx1+brs
      if(ispin.eq.2)vxc=ars*trx2+brs
      uxc=vxc/2.d0
      return
      end
c ************************************************************
      subroutine schrhs(v,k,l,ukl)
      implicit double precision (a-h,o-z)
c integrates the radial schrodinger equation
c on the herman-skillman mesh
      double precision k
      dimension v(561),ukl(561),y(561)
      common/mess/wj(301),dx,nblock,jblock
      common/sc/gr(561),r(561),snlo(561),nbl
      sd(i)=(2.d0*v(i)-k**2+l*(l+1)/r(i)**2)*ukl(i)
      sdd(i)=(2.d0*v(i)-k**2+l*(l+1)/r(i)**2)*y(i)
c initial values have been stored at ukl(1)...ukl(4)
       nb1=nbl+1
      hh=dx
      do 10 i=5,nb1
      ukl(i)=ukl(i-1)+ukl(i-3)-ukl(i-4)
     1+hh**2/4.d0*(5.d0*sd(i-1)+2.d0*sd(i-2)+5.d0*sd(i-3))
  10  ukl(i)=2.d0*ukl(i-1)-ukl(i-2)+
     1hh**2/12.d0*(sd(i)+10.d0*sd(i-1)+sd(i-2))
c set initial values for the next block
      do 20 j=2,jblock
      hh=2.d0*hh
      jj=(j-1)*nbl
      y(jj+1)=ukl(jj+1)
      y(jj)=ukl(jj-1)
      y(jj-1)=ukl(jj-3)
      y(jj-2)=ukl(jj-5)
      do 30 jk=2,nb1
      i=jj+jk
      y(i)=y(i-1)+y(i-3)-y(i-4)+
     1hh**2/4.d0*(5.d0*sdd(i-1)+2.d0*sdd(i-2)+5.d0*sdd(i-3))
      y(i)=2.d0*y(i-1)-y(i-2)+
     1hh**2/12.d0*(sdd(i)+10.d0*sdd(i-1)+sdd(i-2))
  30  ukl(i)=y(i)
  20  continue
      return
      end