subroutine divdif(n,x,func,second)
      implicit none
C 
C Compute numerically the second derivatives of a function of n variables.
C The lower triangle of the derivative matrix is stored in the vector second. 
C
      integer n,nf,i,ii,j,icant,uiparm
      double precision x(n),xi,xj,f,fxp,fxm,fxpym,fxmyp,dwaf,
     & delta /5.0d-4/, ddelta /2.5d-7/
      double precision urparm
      double precision second(n*(n+1)/2)
      external func,ufparm
      call func(n,x,nf,f,uiparm,urparm,ufparm)
      dwaf=2*f
      do i=1,n
        ii=icant(i,i)
c        print*, i,i," ii",ii
        xi=x(i)
        x(i)=xi+delta
        call func(n,x,nf,fxp,uiparm,urparm,ufparm)
        x(i)=xi-delta
        call func(n,x,nf,fxm,uiparm,urparm,ufparm)
        x(i)=xi
        second(ii)=(fxp+fxm-dwaf)/ddelta
      enddo
      do i=1,n
        do j=1,i-1
          xi=x(i)       
          xj=x(j)
          ii=icant(i,j)
c          print *,i,j," ii",ii
          x(i)=xi+delta
          x(j)=xj-delta
          call func(n,x,nf,fxpym,uiparm,urparm,ufparm)
          x(i)=xi-delta
          x(j)=xj+delta
          call func(n,x,nf,fxmyp,uiparm,urparm,ufparm)
          second(ii)=0.5d0*((second(icant(i,i))+second(icant(j,j)))
     &      -(fxpym+fxmyp-dwaf)/ddelta)
          x(i)=xi
          x(j)=xj
        enddo
      enddo
      return
      end
c------------------------------------------------------------------------------
      subroutine func(n,x,nf,f,uiparm,urparm,ufparm)
      implicit none
      integer n,nf,uiparm
      double precision x(2),urparm,f
      external ufparm
c      f = x(1)**2+x(2)**2
      f=dcos(x(1)-x(2))
      return
      end
C-----------------------------------------------------------------------
      INTEGER FUNCTION ICANT(I,J)
      implicit none
      integer i,j
      IF (I.GE.J) THEN
        ICANT=(I*(I-1))/2+J
      ELSE
        ICANT=(J*(J-1))/2+I
      ENDIF
      RETURN
      END
c------------------------------------------------------------------------
      subroutine ufparm
      return
      end