Sub GaussJordan(matrix() As Double,rhs() As Double,ans() As Double)
    Dim As Integer 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 Integer=1 To n 'take copies
        r(c)=rhs(c)
        For d As Integer=1 To n
            b(c,d)=matrix(c,d)
        Next d
    Next c
    #macro pivot(num)
    For p1 As Integer  = num To n - 1
        For p2 As Integer  = p1 + 1 To n
            If Abs(b(p1,num))<Abs(b(p2,num)) Then
                Swap r(p1),r(p2)
                For g As Integer=1 To n
                    Swap b(p1,g),b(p2,g)
                Next g
            End If
        Next p2
    Next p1
    #endmacro
    For k As Integer=1 To n-1
        pivot(k)              'full pivoting
        For row As Integer =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 Integer=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 Integer=n To 1 Step -1
        ans(z)=r(z)/b(z,z)
        For j As Integer = n To z+1 Step -1
            ans(z)=ans(z)-(b(z,j)*ans(j)/b(z,z))
        Next j
        Next    z
    End Sub

    'Interpolate through points.
    Sub Interpolate(x_values() As Double,y_values() As Double,p() As Double)
        Var n=Ubound(x_values)
        Redim p(0):Redim p(1 To n)
        Dim As Double matrix(1 To n,1 To n),rhs(1 To n)
        For a As Integer=1 To n
            rhs(a)=y_values(a)
            For b As Integer=1 To n
                matrix(a,b)=x_values(a)^(b-1)
            Next b
        Next a
        'Solve the linear equations
        GaussJordan(matrix(),rhs(),p())
    End Sub

 '======================== SET UP THE POINTS ===============

    Dim As Double x(1 To ...)={0,1,2,3,4,5,6,7,8,9,10}
    Dim As Double y(1 To ...)={1,6,17,34,57,86,121,162,209,262,321}

    Redim As Double Poly(0)
    'Get the polynomial Poly()
    Interpolate(x(),y(),Poly())

    'print coefficients to console
    print "Polynomial Coefficients:"
    print
    For z As Integer=1 To Ubound(Poly)
        If z=1 Then
            Print "constant term  ";tab(20);Poly(z)
        Else
            Print tab(8); "x^";z-1;" =  ";tab(20);Poly(z)
        End If
    Next z

    sleep