docs: def step(self, tau=1.0, state_fb=None):
          """Run the DMP system for a single timestep.
  
         tau float: scales the timestep
                    increase tau to make the system execute faster
         state_fb np.array: optional system feedback
          """
  
          # run canonical system
          cs_args = {'tau':tau,
                     'error_coupling':1.0}
          if state_fb is not None: 
              # take the 2 norm of the overall error
              state_fb = state_fb.reshape(1,self.dmps)
              dist = np.sqrt(np.sum((state_fb - self.y)**2))
              cs_args['error_coupling'] = 1.0 / (1.0 + 10*dist)
          x = self.cs.step(**cs_args)
  
          # generate basis function activation
          psi = self.gen_psi(x)
  
          for d in range(self.dmps):
  
              # generate the forcing term
              f = self.gen_front_term(x, d) * \
                  (np.dot(psi, self.w[d])) / np.sum(psi) if self.bfs > 0. else 0.
  
              # DMP acceleration
              self.ddy[d] = (self.ay[d] * 
                       (self.by[d] * (self.goal[d] - self.y[d]) - \
                       self.dy[d]/tau) + f) * tau
              self.dy[d] += self.ddy[d] * tau * self.dt * cs_args['error_coupling']
              self.y[d] += self.dy[d] * self.dt * cs_args['error_coupling']
  
          return self.y, self.dy, self.ddy