The paper introduces a method for calculating sensitivity functions of the first and second orders of phase coordinates that do not require the integration of cumbersome chain-related systems of differential equations. The general and particular solutions of the vector differential equation of motion, written in the Cauchy form, are expressed in terms of the fundamental matrix. With the help of formal transformations, the vector of phase coordinates is represented as a convergent integro-power series with respect to the matrix that determines the variations of the elements of the matrix of system equation coefficients. In this work, the system is called the variation matrix. Then the relations obtained are transformed by special operations to an explicit form with respect to these variations up to and including the quadratic approximation. Within the framework of the matrix apparatus, there are similar expansions in terms of variations of the external influence and initial conditions. To exclude numerous "parasitic" operations of multiplication by zero in calculations, it is proposed to use special operations of matrix algebra. The original and inverse fundamental matrices are treated as multiplicative integrals, which ensures their simple calculation in time using recurrent formulas. The resulting apparatus is built entirely on matrix operations, which ensures simple machine implementation and versatility.