def solution(func: Callable[[int], int] = question_function, order: int = 10) -> int:
    data_points: list[int] = [func(x_val) for x_val in range(1, order + 1)]

    polynomials: list[Callable[[int], int]] = [
        interpolate(data_points[:max_coeff]) for max_coeff in range(1, order + 1)
    ]

    ret: int = 0
    poly: Callable[[int], int]
    x_val: int

    for poly in polynomials:
        x_val = 1
        while func(x_val) == poly(x_val):
            x_val += 1

        ret += poly(x_val)

    return ret