I need to solve a number of non linear equations using VBA (specifically this is for estimation of the equilibrium composition of an gas mixture using Gibbs free energy minimization which will have at least 10 non lin. equations to solve). The problem is that the code I have (from "EXCEL FOR SCIENTISTS & ENGINEERS by E.Joseph Billo) using gaussian jordan elimination/matrix pivoting, it produces a zero pivot & thus an error when the matrix elements are divided by the pivot.

Can anyone suggest how I augment the code to overcome this problem (division by zero pivot term) or has anyone else a better code for non linear equations solver. CODE BELOW

''''''''''''''''''''CODE'''''''''''''''''''''''''''''''
Option Explicit
Option Base 1

Const ConvergenceTolerance = 0.00000001
Const IncermentNumericalDifferentiation = 0.000000001
Const Iterations = 50
Const expon = 2.718281828



Function SimultEqNL(equations, Variables, constants)
' Newton iteration method to find roots of nonlinear simultaneous equations
' Example:
'  w^3 + 2w^2 + 3w + 4 =   +12.828
'  w.x + x.y  + y.z    =   -3.919
'  w^2 + 2w.x + x^2    =   +1
'  w   + x    + y  - z =   -3.663
'
'  WHERE:  constants = [12.828,-3.919,1,-3.663];

'On Error Resume Next

Dim i As Integer, j As Integer, k As Integer, N As Integer
Dim Nlterations As Integer
Dim R As Integer, C As Integer
Dim VarAddr() As String, FormulaString() As String
Dim con() As Double, A() As Double, B() As Double
Dim V() As Double
Dim Y1 As Double, Y2 As Double
Dim tolerance As Double, incr As Double

N = equations.Rows.count
k = Variables.Rows.count

If k = 1 Then k = Variables.Columns.count
If k <> N Then SimultEqNL = CVErr(xlErrRef): Exit Function
' Use the CVErr function to create user-defined errors in user-created procedures.
' For example, if you create a function that accepts several arguments and normally returns a string,
' you can have your function evaluate the input arguments to ensure they are within acceptable range.
' If they are not, it is likely your function will not return what you expect.
' In this event, CVErr allows you to return an error number that tells you what action to take.


ReDim VarAddr(N), FormulaString(N), V(N), con(N)
ReDim A(N, N + 1), B(N, N + 1)

tolerance = ConvergenceTolerance         'Convergence criterion.
incr = IncermentNumericalDifferentiation 'Increment for numerical differentiation.
Nlterations = Iterations

For i = 1 To N
VarAddr(i) = Variables(i).Address
' i.e. VarAddr(1) = $A$11
Next i

' Initial values
For i = 1 To N
con(i) = constants(i).Value
' Put the initial guesses into vector V()
V(i) = Variables(i).Value: If V(i) = 0 Then V(i) = 1
Next i

    For j = 1 To Nlterations
            ' Create N x N matrix of partial derivatives.
                For R = 1 To N   ' n = equations.Rows.count
                    For C = 1 To N
                ' Formulastring is formula in which all but one variable in each equation is replaced by current values.
                    FormulaString(R) = Application.ConvertFormula(equations(R).Formula, xlA1, xlA1, xlAbsolute)
                    ' xlA1 = Use xlA1 to return an A1-style reference. xlR1C1 = Use xlR1C1 to return an R1C1-style reference
                    ' xlAbsolute = Convert to absolute row and column style

                    ' ConvertFormula method used to convert cell references from A1 reference style to R1C1 reference style
                            For i = 1 To N
                           ' Debug.Print FormulaString(R)
                            If i <> C Then FormulaString(R) = Application.Substitute(FormulaString(R), VarAddr(i), V(i))
                            ' Substitutes new_text for old_text in a text string.
                            ' FormulaString(R) = the reference to a cell containing text for which you want to substitute characters.
                            ' Replace the address reference with the actual variable value i.e. $B$5^3 + 2*$B$5^2 + 3*$B$5 + 4-$R$5+3-1^3
                            ' i = C means its on the diagonal of the matrix
                        Next i
                    ' V() = vector of current variable values
                    If IsError(Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), V(C) * (1 + incr)))) Then MsgBox "ERROR IS FORMAULA EVALUATION"
                    If IsError(Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), V(C) * (1 - incr)))) Then MsgBox "ERROR IS FORMAULA EVALUATION"
                    
'                    Y2 = Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), (V(C) + incr)))
'                    ' value of the equation at the current variable figure i.e. instead of evaluating the equation at VarAddr(C) evaluate at V(C)*(1+incr)
'                    Y1 = Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), (V(C) - incr)))
'                    A(R, C) = (Y2 - Y1) / (2 * incr)
                    Y2 = Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), V(C) * (1 + incr)))
                    ' value of the equation at the current variable figure i.e. instead of evaluating the equation at VarAddr(C) evaluate at V(C)*(1+incr)
                    Y1 = Evaluate(Application.Substitute(FormulaString(R), VarAddr(C), V(C) * (1 - incr)))
                    A(R, C) = (Y2 - Y1) / (2 * incr * V(C))
                    '
                     ' Derivatives i.e Taylor Series approx (numerical differentiation) central difference:
                     ' y'  =  (y,i+1  -  y,i-1) / 2h  where y' = derivative of y, h = step size
                     ' F(x+h) = F(x)  +  hF'(x)  thus F'(x)  = [ F(x+h)-F(x) ] /  h
                    Next C
                Next R
                
        'Augment matrix of derivatives with vector of constants.
            For R = 1 To N
                FormulaString(R) = Application.ConvertFormula(equations(R).Formula, xlA1, xlA1, xlAbsolute)
                    For C = 1 To N
                        FormulaString(R) = Application.Substitute(FormulaString(R), VarAddr(C), V(C))
                    Next C
                    If IsError(Evaluate(FormulaString(R))) Then MsgBox "ERROR IN FORMULA FO Augment matrix of derivatives"
                A(R, N + 1) = con(R) - Evaluate(FormulaString(R))
            Next R

            For i = 1 To N
                If Abs((A(i, N + 1)) / V(i)) > tolerance Then GoTo Refine
            Next i

            SimultEqNL = Application.Transpose(V)
            Exit Function

Refine:     Call GaussJordan3(N, A, B)
    'Update V values
        For i = 1 To N
        V(i) = V(i) + A(i, N + 1)
        Next i
        'Debug.Print j, "", V(1), V(2), V(3), V(4) ', V(5), V(6), V(7), V(8), V(9)
Next j
' Exit here if no convergence after a specified number of iteration
SimultEqNL = CVErr(xlErrNA)

End Function
'. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Sub GaussJordan3(N, AugMatrix, TempMatrix)
Dim i As Integer, j As Integer, k As Integer, L As Integer, LL As Integer, P As Integer, M As Integer, MM As Integer, MMM As Integer, MMMM As Integer
Dim pivot As Double, temp As Double, arr1() As Variant, arr2() As Variant, x As Integer, y As Integer
Dim determine As Integer, xx As Integer, yy As Integer, cntr As Integer, fudgefactor As Double
    determine = 0
    x = UBound(AugMatrix, 1): xx = UBound(TempMatrix, 1) ' =Number of ROWS
    y = UBound(AugMatrix, 2): yy = UBound(TempMatrix, 2) ' =Number of COLUMNS' =Number of COLUMNS
  '  Debug.Print x, y, xx, yy
    ReDim arr1(1 To 1, 1 To y): ReDim arr2(1 To 1, 1 To y)

For k = 1 To N
        ' Locate largest matrix element, use as pivot.
     pivot = AugMatrix(k, k): P = k
     
     For L = k + 1 To N   ' loop each row
     
        If Abs(AugMatrix(L, k)) < Abs(pivot) Then GoTo EndOfLoop
        
           pivot = AugMatrix(L, k)
        
        P = L
EndOfLoop:     Next L  ' next row
'Debug.Print pivot
        ' Swap rows
        For j = 1 To N + 1
            temp = AugMatrix(k, j)
            AugMatrix(k, j) = AugMatrix(P, j)
            AugMatrix(P, j) = temp
        Next j
        
        ' Normalize pivot row
        For j = 1 To (N + 1)
       'If pivot = 0 Then MsgBox "PIVOT is ZERO"
       If pivot = 0 Then Debug.Print "TITS",
            TempMatrix(k, j) = AugMatrix(k, j) / pivot
        Next j
        
        ' Do the Gauss elimination.
        For i = 1 To N
            If i = k Then GoTo EndOfLoop2
            For j = 1 To N + 1
                TempMatrix(i, j) = AugMatrix(i, j) - AugMatrix(i, k) * TempMatrix(k, j)
            Next j
EndOfLoop2: Next i
        For i = 1 To N
            For j = 1 To N + 1
                AugMatrix(i, j) = TempMatrix(i, j)
            Next j
        Next i
Debug.Print k, pivot
Next k
End Sub