import numpy as np
from numpy.linalg import solve
def polyfit_4(x, f):
    P= np.array([[1, x[0], x[0]**2, x[0]**3],
                 [1, x[1], x[1]**2, x[1]**3],
                 [1, x[2], x[2]**2, x[2]**3],
                 [1, x[3], x[3]**2, x[3]**3]], dtype = float)

    
    A = solve(P,f);
    coeff=np.array([A[0]-A[1],A[1]-A[2],A[2]-A[3],A[3]], dtype=float)
    return coeff