  The solution of a (stochastic) differential equation can be locally
approximated by a (stochastic) expansion. If the vector field of the
differential equation is a polynomial, the corresponding expansion is a linear
combination of iterated integrals of the drivers and can be calculated using
Picard Iterations. However, such expansions grow exponentially fast in their
number of terms, due to their specific algebra, rendering their practical use
limited.
  We present a Mathematica procedure that addresses this issue by
re-parametrising the polynomials and distributing the load in as small as
possible parts that can be processed and manipulated independently, thus
alleviating large memory requirements and being perfectly suited for
parallelized computation. We also present an iterative implementation of the
shuffle product (as opposed to a recursive one, more usually implemented) as
well as a fast way for calculating the expectation of iterated Stratonovich
integrals for Brownian Motion.
