write a Python program which can make Nth order polynomial fit not using numpy

Here's an example Python program to make an Nth order polynomial fit without using numpy:

def polynomial_fit(x, y, degree):
    # Create a Vandermonde matrix of the given degree
    vandermonde = [[xi**i for i in range(degree+1)] for xi in x]
    
    # Solve the linear system to obtain the coefficients of the polynomial
    coefficients = gaussian_elimination(vandermonde, y)
    
    # Create a lambda function to represent the polynomial
    polynomial = lambda xi: sum([coefficients[i]*xi**i for i in range(degree+1)])
    
    return polynomial, coefficients

def gaussian_elimination(A, b):
    # Perform Gaussian elimination with partial pivoting to solve Ax=b
    n = len(A)
    for i in range(n):
        # Find the pivot row and swap with current row
        pivot_row = max(range(i, n), key=lambda j: abs(A[j][i]))
        A[i], A[pivot_row] = A[pivot_row], A[i]
        b[i], b[pivot_row] = b[pivot_row], b[i]
        
        # Eliminate lower rows
        for j in range(i+1, n):
            factor = A[j][i]/A[i][i]
            for k in range(i+1, n):
                A[j][k] -= factor*A[i][k]
            b[j] -= factor*b[i]
    
    # Back-substitute to obtain solution
    x = [0]*n
    for i in range(n-1, -1, -1):
        x[i] = (b[i] - sum([A[i][j]*x[j] for j in range(i+1, n)]))/A[i][i]
    
    return x

# Example usage
x = [1, 2, 3, 4, 5]
y = [2.1, 3.2, 4.3, 5.4, 6.5]
degree = 3

polynomial, coefficients = polynomial_fit(x, y, degree)
print("Polynomial coefficients:", coefficients)
print(f"Polynomial function: {coefficients[0]:.2f} + {coefficients[1]:.2f}x + {coefficients[2]:.2f}x^2 + {coefficients[3]:.2f}x^3")
print("Predicted value at x=6:", polynomial(6))

This program defines two functions: `polynomial_fit` and `gaussian_elimination`. The `polynomial_fit` function takes in three parameters: `x` and `y` are lists representing the input and output data, and `degree` is the degree of the desired polynomial fit. It first creates a Vandermonde matrix of the given degree using the input data, and then solves the linear system `Ax=b` to obtain the coefficients of the polynomial using the `gaussian_elimination` function. Finally, it creates a lambda function representing the polynomial and returns it along with the coefficients.

The `gaussian_elimination` function is a simple implementation of Gaussian elimination with partial pivoting to solve a linear system `Ax=b`. It first performs partial pivoting to find the pivot element of each row, and then eliminates the lower rows using the pivot element. Finally, it performs back-substitution to obtain the solution.

To use the `polynomial_fit` function, simply call it with the input and output data and the desired degree of the polynomial fit. It returns a lambda function representing the polynomial and a list of coefficients. You can then use the lambda function to predict output values at new input values.