I'm an amateur at math: is what is being calculated in the article the same as a Taylor series?
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!
In this case, the author has picked seven points on the sine curve (or its integral or its derivative), and used those seven pieces of information to find a unique polynomial.
If I were to guess, I'd say that if the seven points were on the curve itself (not on its integral or derivative), and they were all very close to some point P, then the resulting polynomial would look a lot like the taylor series expansion around P.
But that's just a guess :)
In fact, with sine modelling, Taylor series beats this approach with lower error. But with this you get to choose your properties so you might get more useful model instead.
Say, if you want to conjoin two sine models got from Taylor series to have it on a full [0, 2pi], you will have derivative discontinuity - a bump - on the joint. But with this model you wouldn't since it's explicitly required to be the inverse number at 0 and at pi/2.