Derivatives of equations required for the Jacobian [vba]

Question

I found the an excel file on https://www.piping-tools.net/Downloads/1.Listado/ named Math. System of nonlinear equations solved with the Newton-Raphson method, in Excel and VBA.xls, I addepted the vba code so not to have to solve the Derivatives of the equation analytically but rather numerical. Part of the added code:

"dF1x = (-F1(x + 2 * h, y) + 8 * F1(x + h, y) - 8 * F1(x - h, y) + F1(x - 2 * h, y)) / 12 / h "

But now i can't seem to remember from where it is that I stole this idea, formula, it looks like something from the Runge–Kutta methods but then only for x and not y et visa versa for the other formula, see below code.

Is this a correct way of calculating the deravative in the intial point (x=5, y=10) and the subsequent following points, and can this formula be extended to be more precise derivative ? And will this method still work for more complex equations?

First function;

Private Function F1(x, y)
    F1 = 4 * x ^ 2 - y ^ 3 + 28         ' = 0
End Function

Second Function;

Private Function F2(x, y)
    F2 = 3 * x ^ 3 + 4 * y ^ 2 - 145 '145     ' = 0
End Function

Math. System of nonlinear equations solved with the Newton-Raphson method, in Excel and VBA.xls;

'Rev. cjc. 12.02.2016
Function Newton_Example_3_1(x, y)                          'Example 3
            
'x = 5
'y = 10
            
On Error GoTo errHandler
'Solution of a system of two equations, using the Newton-Raphson method
'Initial solutions assumed: x and y
'The system equations and four derivatives must be entered in the Input Box below
            
Dim J(0 To 1, 0 To 1), Jinv(0 To 1, 0 To 1) As Single
Dim Old(0 To 1, 0 To 0), F(0 To 1, 0 To 0), JinvF(0 To 1, 0 To 0), Sol(0 To 1, 0 To 0) As Single
Dim Z, Jinvert As Variant

'Iteratio will stop if in two consecutive iterations,
'both x and y values have a change less than the Dstop-value
Dstop = 0.00001
'Iteration counter
N = 0
            
1000:    'Iteration
N = N + 1
           
'____________________Input Box_____________________
'                                                  |
'Matrix  F: System of equations                    |
F(0, 0) = F1(x, y)                 ' = 0           |
F(1, 0) = F2(x, y)                 ' = 0           |
'                                                  |
'Derivatives of both equations                     |
'required for the Jacobian                         |

Const h = 0.0001 ' Step size

'Special added formula
dF1x = (-F1(x + 2 * h, y) + 8 * F1(x + h, y) - 8 * F1(x - h, y) + F1(x - 2 * h, y)) / 12 / h '= 40
dF1y = (-F1(x, y + 2 * h) + 8 * F1(x, y + h) - 8 * F1(x, y - h) + F1(x, y - 2 * h)) / 12 / h

dF2x = (-F2(x + 2 * h, y) + 8 * F2(x + h, y) - 8 * F2(x - h, y) + F2(x - 2 * h, y)) / 12 / h '=225
dF2y = (-F2(x, y + 2 * h) + 8 * F2(x, y + h) - 8 * F2(x, y - h) + F2(x, y - 2 * h)) / 12 / h
'                                                  |
'__________________________________________________|


'Jacobian matrix
J(0, 0) = dF1x
J(0, 1) = dF1y
J(1, 0) = dF2x
J(1, 1) = dF2y

'Inverse matrix of Jacobian
Jinvert = Application.WorksheetFunction.MInverse(J)


'Matrix Old
Old(0, 0) = x
Old(1, 0) = y

'Matrix product Jinv and F
Z = Application.WorksheetFunction.MMult(Jinvert, F)

'Newton-Raphson method equation
Sol(0, 0) = Old(0, 0) - Z(1, 1)
Sol(1, 0) = Old(1, 0) - Z(2, 1)

'New solutions
x = Sol(0, 0)
y = Sol(1, 0)

'Variables change
Dx = Abs(Old(0, 0) - x)
Dy = Abs(Old(1, 0) - y)

If (Dx < Dstop And Dy < Dstop) Then
    GoTo 2000
Else
    GoTo 1000
End If

2000:

Newton_Example_3_1 = Array(x, y)
Exit Function

errHandler:
'Stop
Exit Function
Resume

End Function