How to write a simple program $polyapp(x,y,n)$ that finds the coefficient vector $\mathbf{c}$ for a polynomial interpolant

Question

I am self-learning some numerical algorithms from the book Numerical methods in Scientific computing by Dahlquist/Bjorck. Problem 4.1.2 asks to write a program that finds the coefficient vector $\mathbf{c}$ for a polynomial interpolant. I'd like someone to help verify my solution, as my interpolation polynomial does not properly fit the data.

Question. (a) Write a program $c = polyapp(x,y,n)$ that finds the coefficient vector $c$ for a polynomial in $p \in \mathcal{P}_n$, in a shifted power basis, such that $y_i \approx p(x_i)$ for $i=1:m,m \ge n$, in the least squares sense, or study a program that does almost this.

(b) The following data show the development of the Swedish Gross Domestic Product (GDP), quoted from a table made by a group associated with the Swedish Employer's confederation. (The data are expressed in prices of 1985 and scaled so that the value for 1950 is 100.

x = np.array([1950, 1955, 1960, 1965, 1970, 1975, 1980, 1985, 1990])
y = np.array([100.0, 117.7,139.3,179.3,219.3,249.1,267.5,291.5,326.4])

For the upper pairs of the data, compute and plot $p(x)$, $x \in [1950,2000]$. Mark the given data points. Do this for $m=9$ and for say $n=9$ and then for $n=8:-2:2$. Store the results so that comparisons can be made afterwards.

My Attempt.

Solution. The shifted power basis for $\mathcal{P}_n$ is given by,

$$ \begin{align*} p_1(x) &= 1\\ p_2(x) &= (x - d)\\ p_3(x) &= (x - d)^2\\ \vdots\\ p_n(x) &= (x - d)^n \end{align*} $$

We have $m > n$ data points. The interpolation problem is, $p(x_i) = f(x_i)$ for $i=1:m$. In matrix form, $M_n(\mathbf{p}) \mathbf{c} = f$.

$$ \begin{bmatrix} p_1(x_1) & p_2(x_1) & \ldots & p_n(x_1) \\ p_1(x_2) & p_2(x_2) & \ldots & p_n(x_2) \\ p_1(x_3) & p_2(x_3) & \ldots & p_n(x_3) \\ \vdots\\ p_1(x_m) & p_2(x_1) & \ldots & p_n(x_m) \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ \vdots\\ c_n \end{bmatrix} = \begin{bmatrix} f(x_1)\\ f(x_2)\\ f(x_3)\\ \vdots\\ f(x_m) \end{bmatrix} $$

We are interested to find the coefficient vector $\mathbf{c}$ such that the sum of the squared residuals

$$ S(\mathbf{c}) = \sum_{i=1}^{m} \left[\sum_{j=1}^{n}c_j p_j(x_i) - f(x_i) \right]^2 $$

is minimized. The partial derivative with respect to $c_k$ is,

$$ \begin{align*} \frac{\partial S(\mathbf{c})}{\partial c_k} &= 2\sum_{i=1}^{m}p_k(x_i)\left(\sum_{j=1}^{n}c_j p_j(x_i)-f(x_i)\right)\\ &= 2\sum_{i=1}^{m}p_k(x_i) \left(\sum_{j=1}^{n} c_j p_j(x_i)\right) -2 \sum_{i=1}^{m}p_k(x_i)f(x_i) \end{align*} $$

The function is minimized at the critical point where the gradient is zero.

My polyapp(x,y,n) routine is as follows.

def polyapp(x,y,n):
    m = len(x) # No of data points
    a = min(x) # Determine the interval [a,b]
    b = max(x)
    d = (a + b)/2 
    
    # We assume a shifted power basis for a polynomial p in P_n.
    def p(x,k):
        return (x - d)**k
    
    # We find the coefficient vector c for a polynomial p in P_n, such that 
    # y_i = p(x_i) in a least squares sense. The coefficient vector is given by the solution of a square linear
    # algebraic system.
    
    A = np.zeros((n,n))
    
    # The term multiplying c_j in the k'th equation is sum_{i=1}^m p_k(x_i) p_j(x_i)
    
    for k in range(0,n):
        for j in range(0,n):
            for i in range(0,m):
                A[k,j] = A[k,j] + p(x[i],k) * p(x[i],j)
                
    
    # Construct the right hand side vector 
    b = np.zeros(m)
    
    for k in range(0,n):
        for i in range(0,m):
            b[k] = b[k] + y[i]*p(x[i],k)
        b[k] = b[k] * n
            
    c = linalg.solve(A,b)
    return A,c

def polyval(x,c,d,n):
    y = 0
    for j in range(0,n):
        y = y + c[j]*((x - d)**j)
    return y
```