  We propose a fast method for high order approximations of the solution of the
Cauchy problem for the linear non-stationary Stokes system in $R^3$ in the
unknown velocity $\bf u$ and kinematic pressure $P$. The density ${\bf f}({\bf
x},t)$ and the divergence vector-free initial value ${\bf g }(\bf x)$ are
smooth and rapidly decreasing as $|\bf x|$ tends to infinity. We construct the
vector $\bf u$ in the form $\bf u=\bf u_1+\bf u_2$ where $\bf u_1$ solves a
system of homogeneous heat equations and $\bf u_2$ solves a system of
non-homogeneous heat equations with right-hand side ${\bf f}-\nabla P$.
Moreover $P=-{\cal L}( \nabla\cdot \bf f)$ where $\cal L$ denotes the harmonic
potential. Fast semi-analytic cubature formulas for computing the harmonic
potential and the solution of the heat equation based on the approximation of
the data by functions with analitically known potentials are considered. In
addition, the gradient $\nabla P$ can be approximated by the gradient of the
cubature of $P$, which is a semi-analytic formula too. We derive fast and
accurate high order formulas for the approximation of $\bf u_1,\bf u_2$, $P$
and $\nabla P$. The accuracy of the method and the convergence order $2,4,6$
and $8$ are confirmed by numerical experiments.
