docs: def calculate_A_b(eqs, unfolding, verbose=0):
      u"""Calculate the equations in matrix form.
  
      >>> from sympy import symbols, pprint, I
      >>> rho = define_density_matrix(2, explicitly_hermitian=True,
      ...                             normalized=True)
  
      >>> Omega = symbols("Omega")
      >>> delta = symbols("delta", real=True)
      >>> hbar = symbols("hbar", positive=True)
      >>> H = hbar*Matrix([[0, Omega.conjugate()/2], [Omega/2, -delta]])
  
      >>> Ne = 2
      >>> aux = define_frequencies(Ne, explicitly_antisymmetric=True)
      >>> omega_level, omega, gamma = aux
  
      >>> eqs = I/hbar*(rho*H-H*rho) + lindblad_terms(gamma, rho, 2)
  
      >>> from fast import Unfolding
      >>> unfolding = Unfolding(Ne, True, True, True)
      >>> A, b = calculate_A_b(eqs, unfolding)
      >>> pprint(A, use_unicode=True)
      ⎡ -γ₂₁   im(Ω)  -re(Ω)⎤
      ⎢                     ⎥
      ⎢        -γ₂₁         ⎥
      ⎢-im(Ω)  ─────    -δ  ⎥
      ⎢          2          ⎥
      ⎢                     ⎥
      ⎢               -γ₂₁  ⎥
      ⎢re(Ω)     δ    ───── ⎥
      ⎣                 2   ⎦
  
      >>> pprint(b, use_unicode=True)
      ⎡   0   ⎤
      ⎢       ⎥
      ⎢-im(Ω) ⎥
      ⎢───────⎥
      ⎢   2   ⎥
      ⎢       ⎥
      ⎢ re(Ω) ⎥
      ⎢ ───── ⎥
      ⎣   2   ⎦
  
      """
      Ne = unfolding.Ne
      Nrho = unfolding.Nrho
      lower_triangular = unfolding.lower_triangular
      rho = define_density_matrix(Ne, explicitly_hermitian=lower_triangular,
                                  normalized=unfolding.normalized)
  
      rho_vect = unfolding(rho)
      if unfolding.real:
          ss_comp = {rho[i, j]: re(rho[i, j])+I*im(rho[i, j])
                     for j in range(Ne) for i in range(Ne)}
  
      A = []; b = []
      for mu in range(Nrho):
          s, i, j = unfolding.IJ(mu)
          if verbose > 0: print mu
          eq = part_symbolic(eqs[i, j].subs(ss_comp), s)
          eq_new = 0
          row = []
          for nu in range(Nrho):
              variable = rho_vect[nu]
              coefficient = Derivative(eq, variable).doit()
              row += [coefficient]
              eq_new += coefficient*variable
  
          b += [-(eq-eq_new).expand()]
  
          A += [row]
      A = Matrix(A); b = Matrix(b)
      return A, b