  The first order by time partial differential equations are used as models in
applications such as fluid flow, heat transfer, solid deformation,
electromagnetic waves, and others. In this paper we propose the new numerical
method to solve a class of initial-boundary value problems for the PDEs using
one of the known space discrete numerical schemes and a Taylor series expansion
by time. Normally a second order discretization by space is applied while a
first order by time is satisfactory. Nevertheless, in a number of different
problems, discretization both by temporal and by spatial variables is needed of
highest orders, which complicates numerical solution, in some cases
dramatically. Therefore it is difficult to apply the same numerical methods for
the solution of some PDE arrays if their parameters are varying in a wide range
so that in some of them different numerical schemes by time fit the best for
precise numerical solution. The Taylor series based solution strategy for the
non-stationary PDEs in CFD simulations has been proposed here that attempts to
optimise the computation time and fidelity of the numerical solution.
