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]