В Python как добавить прирастающий элемент полинома?

Строю полином Ньютона на основе заданной простой функции синуса. Реализовал промежуточные расчеты, но остановился на финальной стадии - на получении формулы полинома. Возможно, здесь поможет рекурсия, но это неточно. Вот формула полинома

Формула перебирает значения из приведенной ниже таблицы: идем по столбцу иксов и первой строчке рассчитанных дельт (идем до той дельты, какую степень полинома получаем). Например, если степень равна 2, значит, возьмем в первой строке 2 дельты и по столбцу иксов - значения до 2.512 (9 скобок с разностями иксов и будет в последнем блоке полинома)

В формуле — набор постоянных блоков, где циклом перебираются значения, но у меня загвоздка в элементе (x —x_0)**[n]. Это степень полинома n, которую задает пользователь. Здесь [n] значит, что выражение в скобке разворачивается: Использую библиотеку sympy для символьных вычислений: x в формуле будущего полинома должен остаться иксом (как символ, а не его значение). Как реализовать часть повторяющегося в полиноме блока, которая прирастает новой скобкой степени полинома?

Код:

import numpy as np
from sympy import *
import pandas as pd
from scipy.special import factorial

def func(x):
    return np.sin(x)

def poly(order):
    # строим колонки X и Y:
    x_i_list = [round( (0.1*np.pi*i), 4 ) for i in range(0, 11)]
    y_i_list = []
    for x in x_i_list:
        y_i = round( (func(x)), 4 ) 
        y_i_list.append(y_i)

    # получаем дельты:
    n=order
    if n < len(y_i_list):
        result = [ np.diff(y_i_list, n=d) for d in np.arange(1, len(y_i_list)) ]
        print(result)
    else:
        print(f'Определите порядок полинома меньше {len(y_i_list)}')
        
    # Определяем индекс в столбце иксов исходя из степени полинома:
    delta_index=len(result[order-1])-1
    x_index = delta_index
    h = (x_i_list[x_index] - x_i_list[0]) / n  # считаем h
    
    b=x_i_list[x_index]
    a=x_i_list[0]
    y_0=x_i_list[0]
    string_one = [] # список с дельтами первой строки (включая столбец степени полинома)
    for elen in result:
        string_one.append(round(elen[0], 4))
    
    # создаем список для последующего прохода по иксам 
    x_col_list = []
    for col in  x_i_list:
        if col <= x_i_list[x_index]:
            x_col_list.append(col)
    
    x = Symbol('x') # для символьного представления иксов 
    
    # идем по дельтам первой строки:
    for delta in string_one:
        # идем по столбику иксов 
        for arg in x_col_list:
            for n in range(1, order+1):
                polynom = ( delta/(factorial(n)*h**n) )*(x - arg) # здесь я остановился