The differential equation problem (1) can be solved by introducing finite difference approximations for u'' and u' . In case of quadratic damping one can use a geometric mean to approximate u'|u'| and thereby linearize the equations. The result of using such numerical methods is an algorithm for computing u(t) at discrete points in time. Let u^n be the approximation to u at time t_n=n\Delta t , n=1,2,\ldots , where \Delta t is a (small) time interval. For example, if \Delta t = 0.1 , we find approximations u^1 to u at t=0.1 , u^2 at t=0.2 , u^3 to t=0.3 , and so forth. Any value u^{n+1} can be computed if u^n and u^{n-1} are known (i.e., previously computed). The formula for u^{n+1} is, in case of linear damping f(u')=bu' , \begin{equation} u^{n+1} = \left(2mu^n + (\frac{b}{2}\Delta t - m)u^{n-1} + \Delta t^2(F^n - s(u^n)) \right)(m + \frac{b}{2}\Delta t)^{-1}, \tag{2} \end{equation} where F^n means F(t) evaluated for t=t_n . A special formula must be applied for n=0 : \begin{equation} u^1 = u^0 + \Delta t\, V + \frac{\Delta t^2}{2m}(-bV - s(u^0) + F^0) \thinspace . \tag{3} \end{equation} For quadratic damping we have a slightly different formula, \begin{align} u^{n+1} = & \left( m + b|u^n-u^{n-1}|\right)^{-1}\times\nonumber\\ &\quad \left(2m u^n - mu^{n-1} + bu^n|u^n-u^{n-1}| + \Delta t^2 (F^n - s(u^n)) \right), \tag{4} \end{align} and again a special formula for u^1 : \begin{equation} u^1 = u^0 + \Delta t V + \frac{\Delta t^2}{2m}\left(-bV|V| - s(u^0) + F^0\right) \thinspace . \tag{5} \end{equation}
The implementation of the computational algorithm can make use
of an array u to represent u^n as u[n]. The force F(t_n) is
assumed to be available as an array element F[n]. The following
Python function computes u given an array t
with time points t_0,t_1,\ldots ,
the initial displacement I, mass m, damping parameter b,
restoring force s(u), environmental forces F as an array
(corresponding to t).