// that's how you get the formula for the cubic interpolant
from sympy import *

y1, y2, y3, y4, a, b, c, d = symbols('y1 y2 y3 y4 a b c d')

print(solve([
    d - y1,
    a * (1/3)**3 + b * (1/3)**2 + c * 1/3 + d - y2,
    a * (2/3)**3 + b * (2/3)**2 + c * 2/3 + d - y3,
    a + b + c + d - y4,
    ], (a, b, c, d)))