Comment by EsportToys

Closed-form (non-iterative) PID solution:

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

   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

For anticipation, just add an extra initial velocity in the opposite direction and let the closed-form solution handle the time evolution. The main trick here is to keep both position and velocity as state. There is no need to “step through the simulation”.