How to use the sympy.zeros function in sympy

def _compute_smat(self):
        """
        Compute the transformation matrix in Screw (6x6) matrix form.
        """
        self._smat = Screw6(
            tl=self.rot, tr=tools.skew(self.trans),
            bl=sympy.zeros(3, 3), br=self.rot
        )

def grad_2d(self,tens):
        self.coord(self.symb[0],self.symb[1],self.symb[2])
        retens = Matrix(tens.diff(self.symb[0]))
        ct = 1
        for sym in self.symb[1:]:
            retens = retens.row_insert(ct, tens.diff(sym))
            ct += 1
        retens = self.invA*retens
        retens = simplify(transpose(self.rot.inv())*retens)
        reeye=sp.zeros(3, 3)
        ssx=tens
        reeye[0,1]=-ssx[0,1]/self.symb[0]
        reeye[1,1]=ssx[0,0]/self.symb[0]
        return retens.transpose()+reeye
    def coord(self,x1,x2,x3):

def _D_matrix(self):

        D = sym.zeros(len(self.voltage_sources), len(self.voltage_sources))

        return D

# Determine which branch currents are needed.
        self.unknown_branch_currents = []

        for elt in self.elements.values():
            if elt.need_branch_current:
                self.unknown_branch_currents.append(elt.name)
            if elt.need_extra_branch_current:
                self.unknown_branch_currents.append(elt.name + 'X')

        # Generate stamps.
        num_nodes = len(self.node_list) - 1
        num_branches = len(self.unknown_branch_currents)

        self._G = sym.zeros(num_nodes, num_nodes)
        self._B = sym.zeros(num_nodes, num_branches)
        self._C = sym.zeros(num_branches, num_nodes)
        self._D = sym.zeros(num_branches, num_branches)

        self._Is = sym.zeros(num_nodes, 1)
        self._Es = sym.zeros(num_branches, 1)

        # Iterate over circuit elements and fill in matrices.
        for elt in self.elements.values():
            elt.stamp(self)

        # Augment the admittance matrix to form A matrix.
        self._A = self._G.row_join(self._B).col_join(self._C.row_join(self._D))
        # Augment the known current vector with known voltage vector
        # to form Z vector.
        self._Z = self._Is.col_join(self._Es)

def create_u_operator(u, transpose=False):
    dim = u.shape[0]
    op_u = s.zeros(dim * dim, dim)
    if transpose:
        for ir in range(dim):
            for ic in range(dim):
                op_u[dim*ir+ic,ic] = u[ir]
    else:
        for ii in range(dim):
            op_u[dim*ii:dim*(ii+1),ii] = u
    return op_u

def weight_function_points(self, func_points, direction, order, block, side):
        """ Multiply the function points with their weightings to form the derviatives."""
        if order == 1:
            h = S.One  # The division of delta is now applied in central derivative to reduce the divisions
            if side == 1:
                h = -S.One*h  # Modify the first derivatives for side ==1
            weighted = zeros(4, 1)
            for i in range(4):
                weighted[i] = h*(self.bc4_coefficients[i, :]*func_points[:, i])
        elif order == 2:
            h_sq = S.One  # The division of delta**2 is now applied in central derivative to reduce the divisions
            weighted = zeros(2, 1)
            for i in range(2):
                weighted[i] = h_sq*(self.bc4_2_coefficients[i, :]*func_points[0:5, i])
        else:
            raise NotImplementedError("Only 1st and 2nd derivatives implemented")
        return weighted

:param derivative_fn:
    :return:
    """

    # Handle the (0)-grade case
    if isinstance(f, sympy.Expr):
        n = len(basis)
        df = [0]*len(basis)
        df = sympy.MutableDenseNDimArray(df)
        for ii in range(n):
            df[ii] = derivative_fn(f, basis[ii])

    # Handle the (1+)-grade cases
    if isinstance(f, sympy.Array) or isinstance(f, sympy.NDimArray):
        n = (len(basis),) + f.shape
        df = sympy.MutableDenseNDimArray(sympy.zeros(*n))
        if len(n) == 2:
            for ii in range(df.shape[0]):
                for jj in range(df.shape[1]):
                    if ii == jj:
                        df[ii, jj] = 0

                    if ii < jj:
                        df[ii, jj] += derivative_fn(f[jj], basis[ii])
                        df[ii, jj] += -derivative_fn(f[ii], basis[jj])
                        df[jj, ii] += -derivative_fn(f[jj], basis[ii])
                        df[jj, ii] += derivative_fn(f[ii], basis[jj])

        if len(n) > 2:
            raise NotImplementedError

        # t = [range(d) for d in df.shape]

def create_vector(name, n_ep, dim):
    """ordering is DOF-by-DOF"""
    vec = s.zeros(dim * n_ep, 1)
    for ii in range(dim):
        for ip in range(n_ep):
            vec[n_ep*ii+ip,0] = '%s%d%d' % (name, ii, ip)
    return vec

deg_f, deg_g are the degrees of the divident and of the
    divisor polynomials respectively, deg_g > deg_f.

    The coefficients of the divident poly are the elements
    in row2 and those of the divisor poly are the elements
    in row1.

    col_num defines the number of columns of the matrix M.

    '''
    if deg_g - deg_f >= 1:
        print('Reverse degrees')
        return

    m = zeros(deg_f - deg_g + 2, col_num)

    for i in range(deg_f - deg_g + 1):
        m[i, :] = rotate_r(row1, i)
    m[deg_f - deg_g + 1, :] = row2

    return m

def mass_matrix_full(self):
        """Augments the coefficients of qdots to the mass_matrix."""

        if self.eom is None:
            raise ValueError('Need to compute the equations of motion first')
        n = len(self.q)
        m = len(self.coneqs)
        row1 = eye(n).row_join(zeros(n, n + m))
        row2 = zeros(n, n).row_join(self.mass_matrix)
        if self.coneqs:
            row3 = zeros(m, n).row_join(self._m_cd).row_join(zeros(m, m))
            return row1.col_join(row2).col_join(row3)
        else:
            return row1.col_join(row2)