полиномиальная поверхность подходит numpy

Это по своей сути численно плохо обусловлено, но вы можете сделать что-то вроде этого:

import numpy as np

x = np.random.randn(500)
y = np.random.randn(500)
z = np.random.randn(500) # Dependent variable

v = np.array([np.ones(500), x, y, x**2, x * y, y**2])

coefficients, residues, rank, singval = np.linalg.lstsq(v.T, z)

Чем больше терминов вы добавляете, тем хуже становится численно. Вы уверены, что вам нужен полиномиальный интерполант?

Есть и другие основания для полиномов, для которых матрица значений не так плохо обусловлена, но я не помню, как они называются; однако любой учебник по численному анализу для колледжа должен иметь этот материал.

Вы можете использовать комбинацию polyvander2d и polyval2d, но вам нужно будет выполнить подгонку самостоятельно, используя выходные данные матрицы проектирования из polyvander2d, возможно, включая масштабирование и тому подобное. Должна быть возможность построить класс Polynomial2d из этих инструментов.