A general algorithm for the numerical solution of systems of the type (1) Mẍ\(+C\dot x+Kx=f\) or (2) Mẍ\(+C\dot xP(x)=f\) in a single time interval [0,\(\Delta\) t] is presented. The order of equation (1) or (2) may be one when the matrix \(M=0\). The solution x is approximated by a polynomial \[ y(t)=\sum^{p-1}_{q=0}x^ q_ n\frac{t^ q}{q!}+\alpha_ n^{(p)\quad}\frac{t^ p}{p!} \] where the values \(x_ n\) of the function x and its derivatives at the starting point \(t=0\) are given and the vector \(\alpha_ n^{(p)}\) is determined. The authors obtain different algorithms using different p and different methods to evaluate the \(\alpha_ n^{(p)}\). The presented algorithm is easily programmed and covers most of the currently used schemes as well as presenting many new possibilities. The authors consider also some particular forms of the general algorithm.