      SUBROUTINE SNDOFD(NR,N,XPLS,FCN,FPLS,A,SX,RNOISE,STEPSZ,ANBR)
C PURPOSE
C -------
C FIND SECOND ORDER FORWARD FINITE DIFFERENCE APPROXIMATION "A"
C TO THE SECOND DERIVATIVE (HESSIAN) OF THE FUNCTION DEFINED BY THE SUBP
C "FCN" EVALUATED AT THE NEW ITERATE "XPLS"
C
C FOR OPTIMIZATION USE THIS ROUTINE TO ESTIMATE
C 1) THE SECOND DERIVATIVE (HESSIAN) OF THE OPTIMIZATION FUNCTION
C    IF NO ANALYTICAL USER FUNCTION HAS BEEN SUPPLIED FOR EITHER
C    THE GRADIENT OR THE HESSIAN AND IF THE OPTIMIZATION FUNCTION
C    "FCN" IS INEXPENSIVE TO EVALUATE.
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C XPLS(N)      --> NEW ITERATE:   X[K]
C FCN          --> NAME OF SUBROUTINE TO EVALUATE FUNCTION
C FPLS         --> FUNCTION VALUE AT NEW ITERATE, F(XPLS)
C A(N,N)      <--  FINITE DIFFERENCE APPROXIMATION TO HESSIAN
C                  ONLY LOWER TRIANGULAR MATRIX AND DIAGONAL
C                  ARE RETURNED
C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
C RNOISE       --> RELATIVE NOISE IN FNAME [F(X)]
C STEPSZ(N)    --> WORKSPACE (STEPSIZE IN I-TH COMPONENT DIRECTION)
C ANBR(N)      --> WORKSPACE (NEIGHBOR IN I-TH DIRECTION)
C
C
      DIMENSION XPLS(N)
      DIMENSION SX(N)
      DIMENSION STEPSZ(N),ANBR(N)
      DIMENSION A(NR,1)
C
C FIND I-TH STEPSIZE AND EVALUATE NEIGHBOR IN DIRECTION
C OF I-TH UNIT VECTOR.
C
      OV3 = 1.0/3.0
      DO 10 I=1,N
        STEPSZ(I)=RNOISE**OV3 * MAX(ABS(XPLS(I)),1./SX(I))
        XTMPI=XPLS(I)
        XPLS(I)=XTMPI+STEPSZ(I)
        CALL FCN(N,XPLS,ANBR(I))
        XPLS(I)=XTMPI
   10 CONTINUE
C
C CALCULATE COLUMN I OF A
C
      DO 30 I=1,N
        XTMPI=XPLS(I)
        XPLS(I)=XTMPI+2.0*STEPSZ(I)
        CALL FCN(N,XPLS,FHAT)
        A(I,I)=((FPLS-ANBR(I))+(FHAT-ANBR(I)))/(STEPSZ(I)*STEPSZ(I))
C
C CALCULATE SUB-DIAGONAL ELEMENTS OF COLUMN
        IF(I.EQ.N) GO TO 25
        XPLS(I)=XTMPI+STEPSZ(I)
        IP1=I+1
        DO 20 J=IP1,N
          XTMPJ=XPLS(J)
          XPLS(J)=XTMPJ+STEPSZ(J)
          CALL FCN(N,XPLS,FHAT)
          A(J,I)=((FPLS-ANBR(I))+(FHAT-ANBR(J)))/(STEPSZ(I)*STEPSZ(J))
          XPLS(J)=XTMPJ
   20   CONTINUE
   25   XPLS(I)=XTMPI
   30 CONTINUE
      RETURN
      END
