c****************************************************************
c     computes gradient by high precision numerical differentiation
c*******************************************************************
c copyright = spellucci 
      subroutine gradfn(x,gradx,err)
      parameter (nx=300)
      implicit real*8 (a-h,o-z)
      logical err
      common/count/icf,icgf,iup,inoup,ires,imod
      real*8 x(nx),gradx(nx),fx(6)
      logical test1,test2,test3,test4
      common/param/epsx,epsg,epsf,test1,test2,test3,test4
      common/mach/epsmach,tolmach,deldif
      common/dim/n
c  epsmach=1.1d-16, tolmach=1.d-292,
c  deldif=.25*exp(log(epsmach)*6.d0/7.d0)

c****************************************************
c     high precision numerical differentiation
c     by sixth order extrapolation
c****************************************************
      err=.false.
      do   i=1,n
        z=x(i)
        xincr=deldif*(dabs(x(i))+1.d0)
        x(i)=z-xincr
        call f(x,fx(1),err)
        if ( err ) return
        x(i)=z+xincr
        call f(x,fx(2),err)
        if ( err ) return
      	xincr=xincr+xincr
        d1=xincr
        x(i)=z-xincr
        call f(x,fx(3),err)
        if ( err ) return
        x(i)=z+xincr
        call f(x,fx(4),err)
        if ( err ) return
        xincr=xincr+xincr
        d2=xincr
        x(i)=z-xincr
        call f(x,fx(5),err)
        if ( err ) return
        x(i)=z+xincr
        call f(x,fx(6),err)
        if ( err ) return
        x(i)=z
        d3=xincr+xincr
        sd1=(fx(2)-fx(1))/d1
        sd2=(fx(4)-fx(3))/d2
        sd3=(fx(6)-fx(5))/d3
        sd3=sd2-sd3
        sd2=sd1-sd2
        sd3=sd2-sd3
        gradx(i)=sd1+.4d0*sd2+sd3/45.d0
      enddo
      icgf=icgf+1
      return
      end
