docs: def store(self, t, step):
          """
          Record the state/algeb values at time t to self.vars
          """
          max_cache = int(self.system.tds.config.max_cache)
          if len(self.vars) >= max_cache > 0:
              self.dump()
              self.vars = list()
              self.t = list()
              self.k = list()
              logger.debug(
                  'varout cache cleared at simulation t = {:g}.'.format(
                      self.system.dae.t))
              self._mode = 'a'
  
          var_data = matrix([self.system.dae.x, self.system.dae.y])
  
          # ===== This code block is deprecated =====
          self.t.append(t)
          self.k.append(step)
          self.vars.append(var_data)
          # =========================================
  
          # clear data cache if written to disk
          if self.np_nrows >= max_cache > 0:
              self.dump_np_vars()
              self.np_vars = np.zeros(self._np_block_shape)
              self.np_nrows = 0
              self.np_t = np.zeros((self._np_block_rows,))
              self.np_k = np.zeros((self._np_block_rows,))
              logger.debug(
                  'np_vars cache cleared at simulation t = {:g}.'.format(
                      self.system.dae.t))
              self._mode = 'a'
  
          # initialize before first-time adding data
          if self.np_nrows == 0:
              self.np_ncols = len(var_data)
              self._np_block_shape = (self._np_block_rows, self.np_ncols)
              self.np_vars = np.zeros(self._np_block_shape)
              self.np_t = np.zeros((self._np_block_rows,))
              self.np_k = np.zeros((self._np_block_rows,))
  
          # adding data to the matrix
          # self.np_vars[self.np_nrows, 0] = t
          self.np_t[self.np_nrows] = t
          self.np_k[self.np_nrows] = step
          self.np_vars[self.np_nrows, :] = np.array(var_data).reshape((-1))
          self.np_nrows += 1
  
          # check if matrix extension is needed
          if self.np_nrows >= self.np_vars.shape[0]:
              self.np_vars = np.concatenate([self.np_vars, np.zeros(self._np_block_shape)], axis=0)
              self.np_t = np.concatenate([self.np_t, np.zeros((self._np_block_rows,))], axis=0)
              self.np_k = np.concatenate([self.np_k, np.zeros((self._np_block_rows,))], axis=0)
  
          # remove the post-computed variables from the variable list
          if self.system.tds.config.compute_flows:
              self.system.dae.y = self.system.dae.y[:self.system.dae.m]