Speaking of 2nd-order linear ODEs w/ const factor (a.k.a. mass-spring-damper systems), I wrote a blog post[0] deriving all possible general solutions in a concise matrix that makes it easily implementable in code.

The following is the complete solution in Lua:

   function sprung_response(t,pos,vel,k,c,m)
      local decay = c/2/m
      local omega = math.sqrt(k/m)
      local resid = decay*decay-omega*omega
      local scale = math.sqrt(math.abs(resid))
      local T1,T0 = t , 1
      if resid<0 then
         T1,T0 = math.sin( scale*t)/scale , math.cos( scale*t)
      elseif resid>0 then
         T1,T0 = math.sinh(scale*t)/scale , math.cosh(scale*t)
      end
      local dissipation = math.exp(-decay*t)
      local evolved_pos = dissipation*( pos*(T0+T1*decay) + vel*(   T1      ) )
      local evolved_vel = dissipation*( pos*(-T1*omega^2) + vel*(T0-T1*decay) )
     return evolved_pos , evolved_vel
   end