def as_real_imag(self, deep=True, **hints):
        from sympy.core.symbol import symbols
        from sympy.polys.polytools import poly
        from sympy.core.function import expand_multinomial
        if self.exp.is_Integer:
            exp = self.exp
            re, im = self.base.as_real_imag(deep=deep)
            if re.func == C.re or im.func == C.im:
                return self, S.Zero
            a, b = symbols('a b', cls=Dummy)
            if exp >= 0:
                if re.is_Number and im.is_Number:
                    # We can be more efficient in this case
                    expr = expand_multinomial(self.base**exp)
                    return expr.as_real_imag()

                expr = poly((a + b)**exp) # a = re, b = im; expr = (a + b*I)**exp
            else:
                mag = re**2 + im**2
                re, im = re/mag, -im/mag
                if re.is_Number and im.is_Number:
                    # We can be more efficient in this case
                    expr = expand_multinomial((re + im*S.ImaginaryUnit)**-exp)
                    return expr.as_real_imag()

                expr = poly((a + b)**-exp)

            # Terms with even b powers will be real
            r = [i for i in expr.terms() if not i[0][1] % 2]
            re_part = Add(*[cc*a**aa*b**bb for (aa, bb), cc in r])
            # Terms with odd b powers will be imaginary
            r = [i for i in expr.terms() if i[0][1] % 4 == 1]
            im_part1 = Add(*[cc*a**aa*b**bb for (aa, bb), cc in r])
            r = [i for i in expr.terms() if i[0][1] % 4 == 3]
            im_part3 = Add(*[cc*a**aa*b**bb for (aa, bb), cc in r])

            return (re_part.subs({a: re, b: S.ImaginaryUnit*im}),
            im_part1.subs({a: re, b: im}) + im_part3.subs({a: re, b: -im}))

        elif self.exp.is_Rational:
            # NOTE: This is not totally correct since for x**(p/q) with
            #       x being imaginary there are actually q roots, but
            #       only a single one is returned from here.
            re, im = self.base.as_real_imag(deep=deep)
            if re.func == C.re or im.func == C.im:
                return self, S.Zero
            r = Pow(Pow(re, 2) + Pow(im, 2), S.Half)
            t = C.atan2(im, re)

            rp, tp = Pow(r, self.exp), t*self.exp

            return (rp*C.cos(tp), rp*C.sin(tp))
        else:

            if deep:
                hints['complex'] = False
                return (C.re(self.expand(deep, **hints)),
                        C.im(self.expand(deep, **hints)))
            else:
                return (C.re(self), C.im(self))

 c=symbols('c0:%d'%(Lamda.rows*Lamda.cols))
 C=Matrix(Lamda.rows,Lamda.cols,c)
 #-----------------------------------------
 
 RHS_matrix=Lamda*sI_A +C                            #right hand side of the equation (1)
  
 '''
 -----------Converting equation (1) to a system of linear -----------
 -----------equations, comparing the coefficients of the  -----------
 -----------polynomials in both sides of the equation (1) -----------
 '''
  EQ=[Eq(LHS_matrix[i,j],expand(RHS_matrix[i,j])) for i,j in product(range(LHS_matrix.rows),range(LHS_matrix.cols)) ]
 eq=[]
 
 for equation in EQ:
     RHS=poly((equation).rhs,s).all_coeffs() #simplify necessary ?
     LHS=poly((equation).lhs,s).all_coeffs()
     if len(RHS)>len(LHS):                      # we add zero for each missing coefficient (greater than  the degree of LHS)
         LHS=(len(RHS)-len(LHS))*[0] + LHS
     eq=eq+[Eq(LHS[i],RHS[i]) for i in range(len(RHS))]
 
 SOL=solve(eq,a+c)   # the coefficients of Λ and C
 
 
 #----------substitute the solution in the matrices----------------
 C=C.subs(SOL)     #substitute(SOL)
 Lamda=Lamda.subs(SOL)    #substitute(SOL)
 Js=(Ds+Cs*Ys+Lamda*(B0));
 # for output compatibility
 E=eye(A0.cols)
 if do_test==True: