SymPy can expand series too, one line of code - and you're done. But that would be no fun.

The idea here is: every differentiation or integration over a polynomial results in another polynomial, right? And a polynomial in a point is just a linear equation of its coefficients. You can shove them into a system and make SymPy solve it.

So you can use differential and integral properties of some random function in some points, turn it into a linear system - and have a model of a function as a polynomial. That's it!