I posted a much more concise derivation (and implementation, in just 16 lines) on HN a while back but didn't gain much traction:

https://news.ycombinator.com/item?id=35899215

Here are the 16 lines in full:

   function sprung_response(t,d,v,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*( d*(T0+T1*decay) + v*(   T1      ) )
      local evolved_vel = dissipation*( d*(-T1*omega^2) + v*(T0-T1*decay) )
     return evolved_pos , evolved_vel
   end

And the embedded interactive demo in the blog post:

https://www.desmos.com/calculator/ynbtai98ns

I'm hesitant about reposting as all of my submissions had been inexplicably blocked by HN filters, by the time someone vouches the post it has already been buried.