%load_ext sympy.interactive.ipythonprinting
import numpy as np
from matplotlib import pyplot as plt
from sympy import *
points = [(-1, 0), (0, -2), (0, -1), (1, 1), (2, 0)]
a, b, c, d = symbols("a b c d")
f = lambda x, y: (y - (a*x**3 + b*x**2 + c*x + d))**2
J = sum(f(*p) for p in points)/2
da = diff(J, a)
db = diff(J, b)
dc = diff(J, c)
dd = diff(J, d)
sols = solve([da, db, dc, dd], [a, b, c, d])
a = sols[a].evalf()
b = sols[b].evalf()
c = sols[c].evalf()
d = sols[d].evalf()
xs = np.linspace(-1.5, 2.5, 1000)
ys = a*xs**3 + b*xs**2 + c*xs + d
pxs, pys = zip(*points)
plt.xlim(-1.5, 2.5)
plt.ylim(-2.5, 1.5)
plt.plot(xs, ys)
plt.plot(pxs, pys, ".")