from sympy import solve, Symbol

def getfunc(p1, p2, p3, p4):
    """ Get a point-returning function for a cubic
        curve over four points, with domain [0 - 3].
    """
    # unknowns
    a, b, c, d, e, f, g, h = [Symbol(n) for n in 'abcdefgh']
    
    # solve for coefficients
    xco = solve([a * 0**3 + b * 0**2 + c * 0 + d - p1.x,
                 a * 1**3 + b * 1**2 + c * 1 + d - p2.x,
                 a * 2**3 + b * 2**2 + c * 2 + d - p3.x,
                 a * 3**3 + b * 3**2 + c * 3 + d - p4.x],

                [a, b, c, d])
    
    yco = solve([e * 0**3 + f * 0**2 + g * 0 + h - p1.y,
                 e * 1**3 + f * 1**2 + g * 1 + h - p2.y,
                 e * 2**3 + f * 2**2 + g * 2 + h - p3.y,
                 e * 3**3 + f * 3**2 + g * 3 + h - p4.y],

                [e, f, g, h])
    
    # shorter variable names
    a, b, c, d = [xco[n] for n in (a, b, c, d)]
    e, f, g, h = [yco[n] for n in (e, f, g, h)]
    
    def func(t):
        """ Return a position for given t or velocity (1st derivative) if arg. is True.
        """
        return Point(a * t**3 + b * t**2 + c * t + d,
                     e * t**3 + f * t**2 + g * t + h)

    return func