Many problems in engineering sciences can be described by linear, inhomogeneous, $m$-th order ordinary differential equations (ODEs) with variable coefficients. For this wide class of problems, we here present a new, simple, flexible, and robust solution method, based on piecewise exact integration of local approximation polynomials as well as on averaging local integrals. The method is designed for modern mathematical software providing efficient environments for numerical matrix-vector operation-based calculus. Based on cubic approximation polynomials, the presented method can be expected to perform (i) similar to the Runge-Kutta method, when applied to stiff initial value problems, and (ii) significantly better than the finite difference method, when applied to boundary value problems. Therefore, we use the presented method for the analysis of engineering problems including the oscillation of a modulated torsional spring pendulum, steady-state heat transfer through a cooling web, and the structural analysis of a slender tower based on second-order beam theory. Related convergence studies provide insight into the satisfying characteristics of the proposed solution scheme.