An approach based on matrix formalism for solving differential equations is described. Effective in sense of performance matrix formalism can be tested with less efficient, but accurate traditional algorithm of numerical simulation based on the Runge-Kutta scheme. In both cases the symplectic version of the algorithms are used. The results coincide to analytical calculations, but some disagreements have been identified. The approach implementation is demonstrated in the problem of long-term spin dynamics in electrostatic fields.