def linearize(self, x0: np.array=None, u0: np.array=None) -> List[np.array]:
        """
        Numerical linearization
        :param x0: initial state
        :param u0: initial input
        :return: list of Jacobians
        """
        ss = self.linearize_symbolic()
        ss_eval = []
        ss_subs = {}
        if x0 is None:
            # noinspection PyUnusedLocal
            x0 = self.x.subs(self.x0)[:]
        if u0 is None:
            # noinspection PyUnusedLocal
            u0 = self.u.subs(self.u0)[:]
        # note, we don't substitute y here since
        # all equations should be in terms of x, u
        # if you substitute y, it will resubstitute
        # over x and cause issues
        ss_subs.update({self.x[i]: x0[i] for i in range(len(self.x))})
        ss_subs.update({self.u[i]: u0[i] for i in range(len(self.u))})
        ss_subs.update(self.p0)
        ss_subs.update(self.c0)
        for i in range(len(ss)):
            ss_eval += [np.matrix(ss[i].subs(ss_subs)).astype(float)]
        return ss_eval