code: def linear_ode_to_matrix(eqs, funcs, t, order):
      r"""
      Convert a linear system of ODEs to matrix form
  
      Explanation
      ===========
  
      Express a system of linear ordinary differential equations as a single
      matrix differential equation [1]. For example the system $x' = x + y + 1$
      and $y' = x - y$ can be represented as
  
      .. math:: A_1 X' = A_0 X + b
  
      where $A_1$ and $A_0$ are $2 \times 2$ matrices and $b$, $X$ and $X'$ are
      $2 \times 1$ matrices with $X = [x, y]^T$.
  
      Higher-order systems are represented with additional matrices e.g. a
      second-order system would look like
  
      .. math:: A_2 X'' =  A_1 X' + A_0 X  + b
  
      Examples
      ========
  
      >>> from sympy import Function, Symbol, Matrix, Eq
      >>> from sympy.solvers.ode.systems import linear_ode_to_matrix
      >>> t = Symbol('t')
      >>> x = Function('x')
      >>> y = Function('y')
  
      We can create a system of linear ODEs like
  
      >>> eqs = [
      ...     Eq(x(t).diff(t), x(t) + y(t) + 1),
      ...     Eq(y(t).diff(t), x(t) - y(t)),
      ... ]
      >>> funcs = [x(t), y(t)]
      >>> order = 1 # 1st order system
  
      Now ``linear_ode_to_matrix`` can represent this as a matrix
      differential equation.
  
      >>> (A1, A0), b = linear_ode_to_matrix(eqs, funcs, t, order)
      >>> A1
      Matrix([
      [1, 0],
      [0, 1]])
      >>> A0
      Matrix([
      [1, 1],
      [1,  -1]])
      >>> b
      Matrix([
      [1],
      [0]])
  
      The original equations can be recovered from these matrices:
  
      >>> eqs_mat = Matrix([eq.lhs - eq.rhs for eq in eqs])
      >>> X = Matrix(funcs)
      >>> A1 * X.diff(t) - A0 * X - b == eqs_mat
      True
  
      If the system of equations has a maximum order greater than the
      order of the system specified, a ODEOrderError exception is raised.
  
      >>> eqs = [Eq(x(t).diff(t, 2), x(t).diff(t) + x(t)), Eq(y(t).diff(t), y(t) + x(t))]
      >>> linear_ode_to_matrix(eqs, funcs, t, 1)
      Traceback (most recent call last):
      ...
      ODEOrderError: Cannot represent system in 1-order form
  
      If the system of equations is nonlinear, then ODENonlinearError is
      raised.
  
      >>> eqs = [Eq(x(t).diff(t), x(t) + y(t)), Eq(y(t).diff(t), y(t)**2 + x(t))]
      >>> linear_ode_to_matrix(eqs, funcs, t, 1)
      Traceback (most recent call last):
      ...
      ODENonlinearError: The system of ODEs is nonlinear.
  
      Parameters
      ==========
  
      eqs : list of SymPy expressions or equalities
          The equations as expressions (assumed equal to zero).
      funcs : list of applied functions
          The dependent variables of the system of ODEs.
      t : symbol
          The independent variable.
      order : int
          The order of the system of ODEs.
  
      Returns
      =======
  
      The tuple ``(As, b)`` where ``As`` is a tuple of matrices and ``b`` is the
      the matrix representing the rhs of the matrix equation.
  
      Raises
      ======
  
      ODEOrderError
          When the system of ODEs have an order greater than what was specified
      ODENonlinearError
          When the system of ODEs is nonlinear
  
      See Also
      ========
  
      linear_eq_to_matrix: for systems of linear algebraic equations.
  
      References
      ==========
  
      .. [1] https://en.wikipedia.org/wiki/Matrix_differential_equation
  
      """
      from sympy.solvers.solveset import linear_eq_to_matrix
  
      if any(
          ode_order(eq, func) > order for eq in eqs for func in funcs
      ):
          msg = "Cannot represent system in {}-order form"
          raise ODEOrderError(msg.format(order))
  
      As = []
  
      for o in range(order, -1, -1):
          # Work from the highest derivative down
          syms = [func.diff(t, o) for func in funcs]
  
          # Ai is the matrix for X(t).diff(t, o)
          # eqs is minus the remainder of the equations.
          try:
              Ai, b = linear_eq_to_matrix(eqs, syms)
          except NonlinearError:
              raise ODENonlinearError(
                  "The system of ODEs is nonlinear."
              )
  
          Ai = Ai.applyfunc(expand_mul)
  
          As.append(Ai if o == order else -Ai)
  
          if o:
              eqs = [-eq for eq in b]
          else:
              rhs = b
  
      return As, rhs
  
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 7 for calling a global function with example usage: linear_ode_to_matrix(eqs, funcs, t, order) and returns: As,  rhs