' ******************************************************************
' Numerical Jacobian
' ******************************************************************

SUB NumJacobian (Equations AS SUB(X() AS DOUBLE, F() AS DOUBLE),_
                 X()       AS DOUBLE,_ 
                 Eta       AS DOUBLE,_ 
                 D()       AS DOUBLE)
' ------------------------------------------------------------------
' Compute the jacobian D of a function vector F at point X
' using finite difference approximations.
' Eta is the relative increment used to compute derivatives.
' ------------------------------------------------------------------

  CONST EtaMin = 1.5D-8  ' Min. relative increment (~ MachEp^(1/2))

  DIM AS INTEGER L, N, I, J
  DIM AS DOUBLE  Temp 

  L = LBOUND(X)
  N = UBOUND(X)
  
  DIM AS DOUBLE Delta(L TO N)   ' Increment
  DIM AS DOUBLE Xminus(L TO N)  ' X - Delta
  DIM AS DOUBLE Xplus(L TO N)   ' X + Delta
  DIM AS DOUBLE Fminus(L TO N)  ' F(X - Delta)
  DIM AS DOUBLE Fplus(L TO N)   ' F(X + Delta)

  IF Eta < EtaMin THEN Eta = EtaMin

  FOR I = L TO N
    IF X(I) <> 0 THEN 
      Delta(I) = Eta * ABS(X(I)) 
    ELSE 
      Delta(I) = Eta
    END IF  
    Xplus(I) = X(I) + Delta(I)
    Xminus(I) = X(I) - Delta(I)
  NEXT I

  FOR J = L TO N
    Temp = X(J)
    
    X(J) = Xminus(J)
    Equations(X(), Fminus())
    
    X(J) = Xplus(J)
    Equations(X(), Fplus())
    
    FOR I = L TO N
      D(I, J) = (Fplus(I) - Fminus(I)) / (2 * Delta(J))
    NEXT I
    
    X(J) = Temp
  NEXT J
END SUB

