A method for increasing the order of approximation to an arbitrary natural number by the numerical integration of boundary value problems for inhomogeneous linear ordinary differential equations of various degrees with variable coefficients by the matrix method

The paper includes the well-known matrix method of numerical integration of boundary value problems for inhomogeneous linear ordinary differential equations with variable coefficients, which provides retaining an arbitrary number of Taylor series expansion members of the sought-for solution or, equally, using the Taylor polynomial of arbitrary degree.
The difference boundary value problem approximating the differential boundary value problem is divided into two subtasks: the first subtask includes difference equations, in the construction of which the boundary conditions of the boundary value problem were not used. The second subtask includes difference equations, in the construction of which the boundary conditions of the problem were used.
Based on the earlier results, the method of increasing the order of approximation of the second subtask per unit, and, consequently, of the entire difference boundary problem as a whole is obtained and tested. The earlier findings are as follows:
a) the order of approximation of the first and second subtasks is proportional to the degree of the Taylor polynomial used;
b) the order of approximation of the first subtask depends on the parity or oddness of the degree of the Taylor polynomial used. It turned out that when using the degrees of the Taylor polynomial which are equal to $2m{-}1$ and $2m$, the approximation orders of these two subtasks are the same;
c) the order of approximation of the second subtask coincides with the order of approximation of the first subtask, if the second subtask does not contain the specified values of any derivatives included in the boundary conditions;
d) the presence in the second subtask of at least one derivative value of varying degrees included in the boundary conditions leads to a decrease in the order of approximation per unit in both the second subtask and the entire difference boundary value problem in general.
The theoretical conclusions have been confirmed by numerical experiments.