Sub GaussJordan(matrix() As Double,rhs() As Double,ans() As Double)
    Dim As Long n=Ubound(matrix,1)
    Redim ans(0):Redim ans(1 To n)
    Dim As Double b(1 To n,1 To n),r(1 To n)
    For c As Long=1 To n 'take copies
        r(c)=rhs(c)
        For d As Long=1 To n
            b(c,d)=matrix(c,d)
        Next d
    Next c
    #macro pivot(num)
    For p1 As Long  = num To n - 1
        For p2 As Long  = p1 + 1 To n
            If Abs(b(p1,num))<Abs(b(p2,num)) Then
                Swap r(p1),r(p2)
                For g As Long=1 To n
                    Swap b(p1,g),b(p2,g)
                Next g
            End If
        Next p2
    Next p1
    #endmacro

    For k As Long=1 To n-1
        pivot(k)              'full pivoting
        For row As Long =k To n-1
            If b(row+1,k)=0 Then Exit For
            Var f=b(k,k)/b(row+1,k)
            r(row+1)=r(row+1)*f-r(k)
            For g As Long=1 To n
                b((row+1),g)=b((row+1),g)*f-b(k,g)
            Next g
        Next row
    Next k
    'back substitute
    For z As Long=n To 1 Step -1
        ans(z)=r(z)/b(z,z)
        For j As Long = n To z+1 Step -1
            ans(z)=ans(z)-(b(z,j)*ans(j)/b(z,z))
        Next j
    Next    z
End Sub

dim as double a(1 to 6,1 to 6) = { _
		                  {1.00, 0.00, 0.00,  0.00,  0.00, 0.00}, _
		                  {1.00, 0.63, 0.39,  0.25,  0.16, 0.10}, _
		                  {1.00, 1.26, 1.58,  1.98,  2.49, 3.13}, _
		                  {1.00, 1.88, 3.55,  6.70, 12.62, 23.80}, _
		                  {1.00, 2.51, 6.32, 15.88, 39.90, 100.28}, _
		                  {1.00, 3.14, 9.87, 31.01, 97.41, 306.02} _
	                      }

dim as double b(1 to 6) = { -0.01, 0.61, 0.91, 0.99, 0.60, 0.02 }

redim as double result()
GaussJordan(a(),b(),result())

for n as long=lbound(result) to ubound(result)
    print result(n)
next n
sleep