Bilateral difference method for solving the boundary value problem for an ordinary differential equation

A method is proposed for calculating the bilateral approximations of the solution of the boundary value problem on $[0, 1]$ for the equation $y''+p(x)y'-q(x)y=f(x)$ and the derivative of the solution having the maximum deviation $O(h^2\omega(h)+h^3)$ on $\{kh\}_{k=0}^N$, where $\omega(t)$ is the sum of the continuity moduli of the functions $p''$, $q''$, $f''$, on the set of points $\{kh\}^N_{k=0}$, $h=1/N$ by means of $O(N)$ operations. The data obtained for fairly smooth $p$, $q$, $f$ allow interpolation to be used for calculating the bilateral approximations of the solution and its higher derivatives having the maximum deviation $O(h^3)$ on $[0, 1]$.