Is it from data science course? :) Here is an almost generic solution I did:

%matplotlib inline

import numpy as np;
import math;
import matplotlib.pyplot as plt;

def f(x):
    return np.sin(x / 5) * np.exp(x / 10) + 5 * np.exp(-x / 2)

# approximate at the given points (feel free to experiment: change/add/remove)
points = np.array([1, 4, 10, 15])
n = points.size

# fill A-matrix, each row is 1 or xi^0, xi^1, xi^2, xi^3 .. xi^n
A = np.zeros((n, n))
for index in range(0, n):
    A[index] = np.power(np.full(n, points[index]), np.arange(0, n, 1))

# fill b-matrix, i.e. function value at the given points
b = f(points)

# solve to get approximation polynomial coefficents
solve = np.linalg.solve(A,b)

# define the polynome approximation of the function
def polinom(x): 
    # Yi = solve * Xi where Xi = x^i
    tiles = np.tile(x, (n, 1))
    tiles[0] = np.ones(x.size)
    for index in range(1, n):
        tiles[index] = tiles[index]**index
    return solve.dot(tiles)

# plot the graphs of original function and its approximation
x = np.linspace(1, 15, 100)
plt.plot(x, f(x))
plt.plot(x, polinom(x))

# print out the coefficients of polynome approximating our function
print(solve)

Based on the previous answer, I solved it like this
This is a abstract solution. To obtain an approximation by the necessary degrees, list them in the degrees array.

import numpy as np
import matplotlib.pyplot as plt

degrees = [1, 2]

def f(x):
    return np.sin(x / 5) * np.exp(x / (2 * 5)) + 5 * np.exp(-x / 2)

def main():
    a, b = 1, 15
    num_points = 100
    coeffs_list = []

    X = np.linspace(a, b, num_points)
    y_original = f(X)
    plt.plot(X, y_original, label='f(x)')

    for degree in degrees:
        points = np.linspace(a, b, degree + 1)
        A = np.vander(points, increasing=True)
        b_vec = f(points)

        print(f"Степень многочлена: {degree}")
        print("Матрица системы:")
        print(A)
        print("Вектор правой части:")
        print(b_vec)

        coeffs_list.append(np.linalg.solve(A, b_vec))

    for coeffs in coeffs_list:
        plt.plot(X, np.polyval(coeffs[::-1], X), label=f'g(x), степень {coeffs.size - 1}')

    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend()
    plt.grid()
    plt.show()

if __name__ == '__main__':
    main()