code: def matrix_exp_jordan_form(A, t):
      r"""
      Matrix exponential $\exp(A*t)$ for the matrix *A* and scalar *t*.
  
      Explanation
      ===========
  
      Returns the Jordan form of the $\exp(A*t)$ along with the matrix $P$ such that:
  
      .. math::
          \exp(A*t) = P * expJ * P^{-1}
  
      Examples
      ========
  
      >>> from sympy import Matrix, Symbol
      >>> from sympy.solvers.ode.systems import matrix_exp, matrix_exp_jordan_form
      >>> t = Symbol('t')
  
      We will consider a 2x2 defective matrix. This shows that our method
      works even for defective matrices.
  
      >>> A = Matrix([[1, 1], [0, 1]])
  
      It can be observed that this function gives us the Jordan normal form
      and the required invertible matrix P.
  
      >>> P, expJ = matrix_exp_jordan_form(A, t)
  
      Here, it is shown that P and expJ returned by this function is correct
      as they satisfy the formula: P * expJ * P_inverse = exp(A*t).
  
      >>> P * expJ * P.inv() == matrix_exp(A, t)
      True
  
      Parameters
      ==========
  
      A : Matrix
          The matrix $A$ in the expression $\exp(A*t)$
      t : Symbol
          The independent variable
  
      References
      ==========
  
      .. [1] https://en.wikipedia.org/wiki/Defective_matrix
      .. [2] https://en.wikipedia.org/wiki/Jordan_matrix
      .. [3] https://en.wikipedia.org/wiki/Jordan_normal_form
  
      """
  
      N, M = A.shape
      if N != M:
          raise ValueError(
              "Needed square matrix but got shape (%s, %s)" % (N, M)
          )
      elif A.has(t):
          raise ValueError("Matrix A should not depend on t")
  
      def jordan_chains(A):
          """Chains from Jordan normal form analogous to M.eigenvects().
          Returns a dict with eignevalues as keys like:
              {e1: [[v111,v112,...], [v121, v122,...]], e2:...}
          where vijk is the kth vector in the jth chain for eigenvalue i.
          """
          P, blocks = A.jordan_cells()
          basis = [P[:, i] for i in range(P.shape[1])]
          n = 0
          chains = {}
          for b in blocks:
              eigval = b[0, 0]
              size = b.shape[0]
              if eigval not in chains:
                  chains[eigval] = []
              chains[eigval].append(basis[n : n + size])
              n += size
          return chains
  
      eigenchains = jordan_chains(A)
  
      # Needed for consistency across Python versions
      eigenchains_iter = sorted(
          eigenchains.items(), key=default_sort_key
      )
      isreal = not A.has(I)
  
      blocks = []
      vectors = []
      seen_conjugate = set()
      for e, chains in eigenchains_iter:
          for chain in chains:
              n = len(chain)
              if (
                  isreal
                  and e != e.conjugate()
                  and e.conjugate() in eigenchains
              ):
                  if e in seen_conjugate:
                      continue
                  seen_conjugate.add(e.conjugate())
                  exprt = exp(re(e) * t)
                  imrt = im(e) * t
                  imblock = Matrix(
                      [[cos(imrt), sin(imrt)], [-sin(imrt), cos(imrt)]]
                  )
                  expJblock2 = Matrix(
                      n,
                      n,
                      lambda i, j: imblock
                      * t ** (j - i)
                      / factorial(j - i)
                      if j >= i
                      else zeros(2, 2),
                  )
                  expJblock = Matrix(
                      2 * n,
                      2 * n,
                      lambda i, j: expJblock2[i // 2, j // 2][
                          i % 2, j % 2
                      ],
                  )
  
                  blocks.append(exprt * expJblock)
                  for i in range(n):
                      vectors.append(re(chain[i]))
                      vectors.append(im(chain[i]))
              else:
                  vectors.extend(chain)
                  fun = (
                      lambda i, j: t ** (j - i) / factorial(j - i)
                      if j >= i
                      else 0
                  )
                  expJblock = Matrix(n, n, fun)
                  blocks.append(exp(e * t) * expJblock)
  
      expJ = Matrix.diag(*blocks)
      P = Matrix(N, N, lambda i, j: vectors[j][i])
  
      return P, expJ
  
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 9 for calling a global function with example usage: matrix_exp_jordan_form(A, t) and returns: P,  expJ, chains