A new method for numerical solving of boundary problem for ordinary differential equations with slowly varying coefficients which is aimed at better representation of solutions in the regions of their rapid oscillations or exponential increasing (decreasing) is proposed. It is based on approximation of the solution to find in the form of superposition of certain polynomial- exponential basic functions. The method is studied for the Helmholtz equation in comparison with the standard finite difference method. The numerical tests have shown the convergence of the method proposed. In comparison with the finite difference method the same accuracy is obtained on substantially rarer mesh. This advantage becomes more pronounced, if the solution varies very rapidly.