Polynomial which satisfies integral and two points

Consider two points (x_0, f_0) and (x_1, f_1) 
let p(x) be the degree two polynomial for which
    p(x_0) = f_0
    p(x_1) = f_1
and the integral of p(x) from -1 to 1 is equal to 0

Write a function which accepts two arguments
    1. a length 2 NumPy vector 'x' of floating point values, with 'x[i]' containing the value of x_i,
    2. a length 2 NumPy vector 'f' of floating point values, with 'f[i]' containing the value of f_i,
and which returns
    a length 3 NumPy vector of floating point values containing the power series coefficients, in order from the highest order term to the constant term, for p(x)

I'm not sure where to start. My intial thought would be to have a differential equation P(1)=P(-1) with initial values p(x_0) = f_0 and p(x_1) = f_1, but I'm also having issues with the implementation.

Using sympy, Python's symbolic math library, the problem can be formulated as follows:

from sympy import symbols Eq, solve, integrate

def give_coeff(x, f):
    a, b, c, X = symbols('a, b, c, X')
    F = a * X * X + b * X + c  # we have a second order polynomial
    sol = solve([Eq(integrate(F, (X, -1, 1)), 0),  # the integral should be zero (2/3*a + 2*c)
                 Eq(F.subs(X, x[0]), f[0]),        # filling in x[0] should give f[0]
                 Eq(F.subs(X, x[1]), f[1])],       # filling in x[1] should give f[1]
                (a, b, c))   # solve for a, b and c
    return sol[a].evalf(), sol[b].evalf(), sol[c].evalf()

import numpy as np

coeff = give_coeff(np.array([1, 2]), np.array([3, 4]))
print(coeff)

The code can even be expanded to polynomials of any degree:

from sympy import Eq, solve, symbols, integrate


def give_coeff(x, f):
    assert len(x) == len(f), "x and f need to have the same length"
    degree = len(x)
    X = symbols('X')
    a = [symbols(f'a_{i}') for i in range(degree + 1)]
    F = 0
    for ai in a[::-1]:
        F = F * X + ai
    sol = solve([Eq(integrate(F, (X, -1, 1)), 0)] +
                [Eq(F.subs(X, xi), fi) for xi, fi in zip(x, f)],
                (*a,))
    # print(sol)
    # print(F.subs(sol).expand())
    return [sol[ai].evalf() for ai in a[::-1]]

import numpy as np

coeff = give_coeff(np.array([1, 2]), np.array([3, 4]))
print(coeff)
print(give_coeff(np.array([1, 2, 3, 4, 5]), np.array([3, 4, 6, 9, 1])))

PS: To solve the second degree equation only using numpy, np.linalg.solve can be used to solve the linear system of 3 unknowns with 3 equations. The equations need to be "hand calculated" which is are more error prone and more elaborated to extend to higher degrees.

import numpy as np

def np_give_coeff(x, f):
    # general equation: F = a*X**2 + b*X + c
    # 3 equations:
    #     integral (F, (X, -1, 1)) == 0 or (2/3*a + 2*c) == 0
    #     a*x[0]**2 + b*x[0] + c == f[0]
    #     a*x[1]**2 + b*x[1] + c == f[1]
    A = np.array([[2/3, 0, 2],
                  [x[0]**2, x[0], 1],
                  [x[1]**2, x[1], 1]])
    B = np.array([0, f[0], f[1]])
    return np.linalg.solve(A, B)

coeff = np_give_coeff(np.array([1, 2]), np.array([3, 4]))
print(coeff)