I feel like publishing some previous works from my PhD thesis. I was using Mathematica to build a system of 2N partial differential equations for 2N functions by symbolic spatial Taylor expansion, then numerically integrate them with NDSolve with respect to time.

Mathematica had a strange behaviour for some N orders, stopping to some apparent numerical singularities I was not able to manage at that time. Furthermore, I have an access to a supercomputer, which doesn't offer Mathematica. I would therefore like to rebuild my model in Python.

For the symbolic part, I guess SymPy will handle easily the Taylor expansion part (notice that after order 2 or 3, it took far more than a full page to Mathematica to write the PDE system, should I ask it to). But how do you numerically (Runge-Kutta or whatever) solve the objects coming out of SymPy?