var("x0 x1 p0 p1 dp0 dp1")
var("a0 a1 a2 a3")
p(x) = a0 + a1*x + a2*x^2 + a3*x^3
dp(x) = derivative(p(x), x)

eq0 = (p(x = x0) == p0)
eq1 = (p(x = x1) == p1)
eq2 = (dp(x = x0) == dp0)
eq3 = (dp(x = x1) == dp1)

res = solve([eq0, eq1, eq2, eq3], [a0, a1, a2, a3])

coeff0 = res[0][0].rhs()
coeff1 = res[0][1].rhs()
coeff2 = res[0][2].rhs()
coeff3 = res[0][3].rhs()

def cubInterpol(xx0, xx1, y0, y1, dy0, dy1):
    return (
         coeff0(x0 = xx0, x1 = xx1, p0 = y0, p1 = y1,
                dp0 = dy0, dp1 = dy1) +
         coeff1(x0 = xx0, x1 = xx1, p0 = y0, p1 = y1,
                dp0 = dy0, dp1 = dy1)*x +
         coeff2(x0 = xx0, x1 = xx1, p0 = y0, p1 = y1,
                dp0 = dy0, dp1 = dy1)*x^2 +
         coeff3(x0 = xx0, x1 = xx1, p0 = y0, p1 = y1,
                dp0 = dy0, dp1 = dy1)*x^3
    )

def plotCubInterpol(
    xx0, xx1, y0, y1, dy0, dy1,
    xmin = -8, xmax = 8, ymin = -6, ymax = 6
):
    cubPol = cubInterpol(xx0, xx1, y0, y1, dy0, dy1)
    return plot(
        cubPol, x, 
        xmin = xmin, xmax = xmax, 
        ymin = ymin, ymax = ymax, aspect_ratio = 1,
        color = "blue"
    )             
     
def drawCubicInterpolation(
    xx0, xx1, y0, y1, dy0, dy1,
    xmin = -8, xmax = 8, ymin = -6, ymax = 6
):
    f = plotCubInterpol(
        xx0, xx1, y0, y1, dy0, dy1,
        xmin, xmax, ymin, ymax
    );
    
    nodes = (
        point([xx0, y0], color="red", size=40) + 
        point([xx1, y1], color="red", size=40)
    )
    
    tangents = (
        line([(xx0 - 1, y0 - dy0), (xx0 + 1, y0 + dy0)], color="green") +        
        line([(xx1 - 1, y1 - dy1), (xx1 + 1, y1 + dy1)], color="green")
    )
    return f + nodes + tangents