We propose a fast method for high order approximations of the solution to the Cauchy problem for the linear nonstationary Stokes system in ℝ3 in the unknown velocity u and kinematic pressure P. The density f (x, t) and the divergence-free vector initial value g(x) are smooth and rapidly decreasing as |x| → ∞. We construct the vector u = u1+u2 where u1 solves a system of homogeneous heat equations and u2 solves a system of nonhomogeneous heat equations with right-hand side f −∇P, where P = −[InlineMediaObject not available: see fulltext.] (∇ · f) and [InlineMediaObject not available: see fulltext.] denotes the harmonic potential. Fast semianalytic cubature formulas for computing the harmonic potential and the solution to the heat equation based on the approximation of the data by functions with analitically known potentials are considered. The gradient ∇P can be approximated by the gradient of the cubature of P, which is a semianalytic formula. We derive fast and accurate high order formulas for approximation of u1, u2, P, and ∇P. The accuracy of the method and the convergence order 2, 4, 6, 8 are confirmed by numerical experiments. © 2019, Springer Science+Business Media, LLC, part of Springer Nature.