def mass_spring(x, u, dt=0.05, m=1, c=2, k=1): 
    """Mass-spring-damper system
    x  :: state variables (x1, x2) where
         x1 = position (m)
         x2 = velocity (m/s)
    u  :: input force (N)
    dt :: timestep (s)
    m  :: mass (kg)
    c  :: damping coefficient (Ns/m)
    k  :: spring constant (N/m)
    """
    dxdt = (x[1], (-k*x[0] - c*x[1] + u) / m)
    return x[0] + dt*dxdt[0], x[1] + dt*dxdt[1]