code: def linodesolve(
      A, t, b=None, B=None, type="auto", doit=False, tau=None
  ):
      r"""
      System of n equations linear first-order differential equations
  
      Explanation
      ===========
  
      This solver solves the system of ODEs of the following form:
  
      .. math::
          X'(t) = A(t) X(t) +  b(t)
  
      Here, $A(t)$ is the coefficient matrix, $X(t)$ is the vector of n independent variables,
      $b(t)$ is the non-homogeneous term and $X'(t)$ is the derivative of $X(t)$
  
      Depending on the properties of $A(t)$ and $b(t)$, this solver evaluates the solution
      differently.
  
      When $A(t)$ is constant coefficient matrix and $b(t)$ is zero vector i.e. system is homogeneous,
      the system is "type1". The solution is:
  
      .. math::
          X(t) = \exp(A t) C
  
      Here, $C$ is a vector of constants and $A$ is the constant coefficient matrix.
  
      When $A(t)$ is constant coefficient matrix and $b(t)$ is non-zero i.e. system is non-homogeneous,
      the system is "type2". The solution is:
  
      .. math::
          X(t) = e^{A t} ( \int e^{- A t} b \,dt + C)
  
      When $A(t)$ is coefficient matrix such that its commutative with its antiderivative $B(t)$ and
      $b(t)$ is a zero vector i.e. system is homogeneous, the system is "type3". The solution is:
  
      .. math::
          X(t) = \exp(B(t)) C
  
      When $A(t)$ is commutative with its antiderivative $B(t)$ and $b(t)$ is non-zero i.e. system is
      non-homogeneous, the system is "type4". The solution is:
  
      .. math::
          X(t) =  e^{B(t)} ( \int e^{-B(t)} b(t) \,dt + C)
  
      When $A(t)$ is a coefficient matrix such that it can be factorized into a scalar and a constant
      coefficient matrix:
  
      .. math::
          A(t) = f(t) * A
  
      Where $f(t)$ is a scalar expression in the independent variable $t$ and $A$ is a constant matrix,
      then we can do the following substitutions:
  
      .. math::
          tau = \int f(t) dt, X(t) = Y(tau), b(t) = b(f^{-1}(tau))
  
      Here, the substitution for the non-homogeneous term is done only when its non-zero.
      Using these substitutions, our original system becomes:
  
      .. math::
          Y'(tau) = A * Y(tau) + b(tau)/f(tau)
  
      The above system can be easily solved using the solution for "type1" or "type2" depending
      on the homogeneity of the system. After we get the solution for $Y(tau)$, we substitute the
      solution for $tau$ as $t$ to get back $X(t)$
  
      .. math::
          X(t) = Y(tau)
  
      Systems of "type5" and "type6" have a commutative antiderivative but we use this solution
      because its faster to compute.
  
      The final solution is the general solution for all the four equations since a constant coefficient
      matrix is always commutative with its antidervative.
  
      An additional feature of this function is, if someone wants to substitute for value of the independent
      variable, they can pass the substitution `tau` and the solution will have the independent variable
      substituted with the passed expression(`tau`).
  
      Parameters
      ==========
  
      A : Matrix
          Coefficient matrix of the system of linear first order ODEs.
      t : Symbol
          Independent variable in the system of ODEs.
      b : Matrix or None
          Non-homogeneous term in the system of ODEs. If None is passed,
          a homogeneous system of ODEs is assumed.
      B : Matrix or None
          Antiderivative of the coefficient matrix. If the antiderivative
          is not passed and the solution requires the term, then the solver
          would compute it internally.
      type : String
          Type of the system of ODEs passed. Depending on the type, the
          solution is evaluated. The type values allowed and the corresponding
          system it solves are: "type1" for constant coefficient homogeneous
          "type2" for constant coefficient non-homogeneous, "type3" for non-constant
          coefficient homogeneous, "type4" for non-constant coefficient non-homogeneous,
          "type5" and "type6" for non-constant coefficient homogeneous and non-homogeneous
          systems respectively where the coefficient matrix can be factorized to a constant
          coefficient matrix.
          The default value is "auto" which will let the solver decide the correct type of
          the system passed.
      doit : Boolean
          Evaluate the solution if True, default value is False
      tau: Expression
          Used to substitute for the value of `t` after we get the solution of the system.
  
      Examples
      ========
  
      To solve the system of ODEs using this function directly, several things must be
      done in the right order. Wrong inputs to the function will lead to incorrect results.
  
      >>> from sympy import symbols, Function, Eq
      >>> from sympy.solvers.ode.systems import canonical_odes, linear_ode_to_matrix, linodesolve, linodesolve_type
      >>> from sympy.solvers.ode.subscheck import checkodesol
      >>> f, g = symbols("f, g", cls=Function)
      >>> x, a = symbols("x, a")
      >>> funcs = [f(x), g(x)]
      >>> eqs = [Eq(f(x).diff(x) - f(x), a*g(x) + 1), Eq(g(x).diff(x) + g(x), a*f(x))]
  
      Here, it is important to note that before we derive the coefficient matrix, it is
      important to get the system of ODEs into the desired form. For that we will use
      :obj:`sympy.solvers.ode.systems.canonical_odes()`.
  
      >>> eqs = canonical_odes(eqs, funcs, x)
      >>> eqs
      [[Eq(Derivative(f(x), x), a*g(x) + f(x) + 1), Eq(Derivative(g(x), x), a*f(x) - g(x))]]
  
      Now, we will use :obj:`sympy.solvers.ode.systems.linear_ode_to_matrix()` to get the coefficient matrix and the
      non-homogeneous term if it is there.
  
      >>> eqs = eqs[0]
      >>> (A1, A0), b = linear_ode_to_matrix(eqs, funcs, x, 1)
      >>> A = A0
  
      We have the coefficient matrices and the non-homogeneous term ready. Now, we can use
      :obj:`sympy.solvers.ode.systems.linodesolve_type()` to get the information for the system of ODEs
      to finally pass it to the solver.
  
      >>> system_info = linodesolve_type(A, x, b=b)
      >>> sol_vector = linodesolve(A, x, b=b, B=system_info['antiderivative'], type=system_info['type_of_equation'])
  
      Now, we can prove if the solution is correct or not by using :obj:`sympy.solvers.ode.checkodesol()`
  
      >>> sol = [Eq(f, s) for f, s in zip(funcs, sol_vector)]
      >>> checkodesol(eqs, sol)
      (True, [0, 0])
  
      We can also use the doit method to evaluate the solutions passed by the function.
  
      >>> sol_vector_evaluated = linodesolve(A, x, b=b, type="type2", doit=True)
  
      Now, we will look at a system of ODEs which is non-constant.
  
      >>> eqs = [Eq(f(x).diff(x), f(x) + x*g(x)), Eq(g(x).diff(x), -x*f(x) + g(x))]
  
      The system defined above is already in the desired form, so we do not have to convert it.
  
      >>> (A1, A0), b = linear_ode_to_matrix(eqs, funcs, x, 1)
      >>> A = A0
  
      A user can also pass the commutative antiderivative required for type3 and type4 system of ODEs.
      Passing an incorrect one will lead to incorrect results. If the coefficient matrix is not commutative
      with its antiderivative, then :obj:`sympy.solvers.ode.systems.linodesolve_type()` raises a NotImplementedError.
      If it does have a commutative antiderivative, then the function just returns the information about the system.
  
      >>> system_info = linodesolve_type(A, x, b=b)
  
      Now, we can pass the antiderivative as an argument to get the solution. If the system information is not
      passed, then the solver will compute the required arguments internally.
  
      >>> sol_vector = linodesolve(A, x, b=b)
  
      Once again, we can verify the solution obtained.
  
      >>> sol = [Eq(f, s) for f, s in zip(funcs, sol_vector)]
      >>> checkodesol(eqs, sol)
      (True, [0, 0])
  
      Returns
      =======
  
      List
  
      Raises
      ======
  
      ValueError
          This error is raised when the coefficient matrix, non-homogeneous term
          or the antiderivative, if passed, are not a matrix or
          do not have correct dimensions
      NonSquareMatrixError
          When the coefficient matrix or its antiderivative, if passed is not a
          square matrix
      NotImplementedError
          If the coefficient matrix does not have a commutative antiderivative
  
      See Also
      ========
  
      linear_ode_to_matrix: Coefficient matrix computation function
      canonical_odes: System of ODEs representation change
      linodesolve_type: Getting information about systems of ODEs to pass in this solver
  
      """
  
      if not isinstance(A, MatrixBase):
          raise ValueError(
              filldedent(
                  """\
              The coefficients of the system of ODEs should be of type Matrix
          """
              )
          )
  
      if not A.is_square:
          raise NonSquareMatrixError(
              filldedent(
                  """\
              The coefficient matrix must be a square
          """
              )
          )
  
      if b is not None:
          if not isinstance(b, MatrixBase):
              raise ValueError(
                  filldedent(
                      """\
                  The non-homogeneous terms of the system of ODEs should be of type Matrix
              """
                  )
              )
  
          if A.rows != b.rows:
              raise ValueError(
                  filldedent(
                      """\
                  The system of ODEs should have the same number of non-homogeneous terms and the number of
                  equations
              """
                  )
              )
  
      if B is not None:
          if not isinstance(B, MatrixBase):
              raise ValueError(
                  filldedent(
                      """\
                  The antiderivative of coefficients of the system of ODEs should be of type Matrix
              """
                  )
              )
  
          if not B.is_square:
              raise NonSquareMatrixError(
                  filldedent(
                      """\
                  The antiderivative of the coefficient matrix must be a square
              """
                  )
              )
  
          if A.rows != B.rows:
              raise ValueError(
                  filldedent(
                      """\
                          The coefficient matrix and its antiderivative should have same dimensions
                      """
                  )
              )
  
      if (
          not any(type == "type{}".format(i) for i in range(1, 7))
          and not type == "auto"
      ):
          raise ValueError(
              filldedent(
                  """\
                      The input type should be a valid one
                  """
              )
          )
  
      n = A.rows
  
      # constants = numbered_symbols(prefix='C', cls=Dummy, start=const_idx+1)
      Cvect = Matrix([Dummy() for _ in range(n)])
  
      if b is None and any(
          type == typ for typ in ["type2", "type4", "type6"]
      ):
          b = zeros(n, 1)
  
      is_transformed = tau is not None
      passed_type = type
  
      if type == "auto":
          system_info = linodesolve_type(A, t, b=b)
          type = system_info["type_of_equation"]
          B = system_info["antiderivative"]
  
      if type in ("type5", "type6"):
          is_transformed = True
          if passed_type != "auto":
              if tau is None:
                  system_info = _first_order_type5_6_subs(A, t, b=b)
                  if not system_info:
                      raise ValueError(
                          filldedent(
                              """
                          The system passed isn't {}.
                      """.format(
                                  type
                              )
                          )
                      )
  
                  tau = system_info["tau"]
                  t = system_info["t_"]
                  A = system_info["A"]
                  b = system_info["b"]
  
      intx_wrtt = lambda x: Integral(x, t) if x else 0
      if type in ("type1", "type2", "type5", "type6"):
          P, J = matrix_exp_jordan_form(A, t)
          P = simplify(P)
  
          if type in ("type1", "type5"):
              sol_vector = P * (J * Cvect)
          else:
              Jinv = J.subs(t, -t)
              sol_vector = (
                  P
                  * J
                  * ((Jinv * P.inv() * b).applyfunc(intx_wrtt) + Cvect)
              )
      else:
          if B is None:
              B, _ = _is_commutative_anti_derivative(A, t)
  
          if type == "type3":
              sol_vector = B.exp() * Cvect
          else:
              sol_vector = B.exp() * (
                  ((-B).exp() * b).applyfunc(intx_wrtt) + Cvect
              )
  
      if is_transformed:
          sol_vector = sol_vector.subs(t, tau)
  
      gens = sol_vector.atoms(exp)
  
      if type != "type1":
          sol_vector = [expand_mul(s) for s in sol_vector]
  
      sol_vector = [
          collect(s, ordered(gens), exact=True) for s in sol_vector
      ]
  
      if doit:
          sol_vector = [s.doit() for s in sol_vector]
  
      return sol_vector
  
imports: from sympy.core import Add, Mul, S
  from sympy.core.containers import Tuple
  from sympy.core.exprtools import factor_terms
  from sympy.core.numbers import I
  from sympy.core.relational import Eq, Equality
  from sympy.core.sorting import default_sort_key, ordered
  from sympy.core.symbol import Dummy, Symbol
  from sympy.core.function import expand_mul, expand, Derivative, AppliedUndef, Function, Subs
  from sympy.functions import exp, im, cos, sin, re, Piecewise, piecewise_fold, sqrt, log
  from sympy.functions.combinatorial.factorials import factorial
  from sympy.matrices import zeros, Matrix, NonSquareMatrixError, MatrixBase, eye
  from sympy.polys import Poly, together
  from sympy.simplify import collect, radsimp, signsimp
  from sympy.simplify.powsimp import powdenest, powsimp
  from sympy.simplify.ratsimp import ratsimp
  from sympy.simplify.simplify import simplify
  from sympy.sets.sets import FiniteSet
  from sympy.solvers.deutils import ode_order
  from sympy.solvers.solveset import NonlinearError, solveset
  from sympy.utilities.iterables import connected_components, iterable, strongly_connected_components
  from sympy.utilities.misc import filldedent
  from sympy.integrals.integrals import Integral, integrate
label_desc: Use image node_id 10 for calling a global function with example usage: linodesolve(A, t, b, B, type, doit, tau) and returns: sol_vector