from sympy import *
import numpy as np
h = symbols('h')
x = symbols('x')
xi = symbols('x_i')
# (9.21)
M = Matrix([
[1, xi, xi**2, xi**3],
[0, 1, 2*xi, 3*(xi**2)],
[1, xi + h, (xi + h)**2, (xi + h)**3],
[0, 1, 2*(xi + h), 3*(xi + h)**2]])
M
M_inv = simplify(M.inv())
M_inv
a = simplify(M_inv.T * Matrix([1, x, x**2, x**3]))
a
M_inv.subs(xi, 0)
a.subs(xi, 0)
# TBD