Вариант решения вашего вопроса в одномерном случае:
>>>import sympy
>>>x = sympy.symbols('x')
>>>sympy.interpolate([(-1, 2), (1, 2), (2, 5)], x)
x**2 + 1
А вот в случае нескольких переменных я не знаю как решить подобную задачу "в лоб".
В маткаде есть функции linfit, polyfit, но это уже в относится к аппроксимации.
