diff --git a/doc/changelog.rst b/doc/changelog.rst
index d02ede32..f9ffa796 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -13,13 +13,18 @@ Changelog
 - :bug:`590` Worked around a performance regression in SymPy 1.13 that
   caused ``examples/ipython/LaTeX.ipynb`` (``check('curvi_linear_latex')``) to
   time out after 600 s on SymPy ≥ 1.13.  SymPy PR #26390 added an O(N·M)
-  ``.replace()`` traversal inside ``TR3``/``futrig`` that is a no-op for
-  galgebra's symbolic trig arguments but dominated each of the ~70
-  ``Simp.apply`` calls during ``Ga.build(norm=True)`` for curvilinear
-  coordinates.  The fix uses ``trigsimp(method='old')`` via ``Simp.profile``
-  for the affected example, cutting run time from > 600 s to < 6 s.
+  ``.replace()`` traversal inside ``TR3``/``futrig`` that made simplification
+  of the prolate-spheroidal output stall during display.  The fix uses
+  ``trigsimp(method='old')`` via ``Simp.profile`` for the affected example,
+  cutting its run time from > 600 s to < 6 s.
   A notebook note documents the two cosmetic output differences from the
-  pre-1.13 form; a proper upstream fix is tracked in :issue:`576`.
+  pre-1.13 form.
+
+- :bug:`598` Multivector string and LaTeX display now avoid the same SymPy
+  regression outside that example.  Large expressions with trigonometric and
+  hyperbolic functions under non-integral powers use the bounded
+  ``trigsimp(method='old')`` path.  Algebraic simplification and explicit
+  ``Simp.profile`` modes remain unchanged.
 
 - :support:`589` Added Step 0 to the release-process runbook
   (``doc/dev/release-process.md``): open a release issue before preparing the
diff --git a/examples/LaTeX/curvi_linear_latex.py b/examples/LaTeX/curvi_linear_latex.py
index 77c54de2..7be83f5a 100644
--- a/examples/LaTeX/curvi_linear_latex.py
+++ b/examples/LaTeX/curvi_linear_latex.py
@@ -189,7 +189,7 @@ def main():
     from sympy import trigsimp
     from galgebra.metric import Simp
 
-    orig_modes = Simp.modes[:]
+    orig_modes = Simp.modes
     Simp.profile([lambda e: trigsimp(e, method='old')])
     try:
         derivatives_in_spherical_coordinates()
diff --git a/galgebra/_utils/simplify.py b/galgebra/_utils/simplify.py
new file mode 100644
index 00000000..921027fd
--- /dev/null
+++ b/galgebra/_utils/simplify.py
@@ -0,0 +1,76 @@
+"""Compatibility helpers for simplification across SymPy releases."""
+
+import re
+
+import sympy
+from sympy import preorder_traversal, simplify, trigsimp
+from sympy.functions.elementary.hyperbolic import HyperbolicFunction
+from sympy.functions.elementary.trigonometric import TrigonometricFunction
+
+
+def _major_minor(version):
+    """Return the leading major and minor numbers from a version string."""
+    match = re.match(r'^(\d+)\.(\d+)', version)
+    if match is None:
+        return (0, 0)
+    return tuple(map(int, match.groups()))
+
+
+_SYMPY_MAJOR_MINOR = _major_minor(sympy.__version__)
+
+# SymPy 1.13's gh-26390 added a nested replace traversal to the FU
+# simplifier. Match only the observed two-term prolate radical shape. A
+# numerical tree-cost heuristic admitted benign expressions whose unrelated
+# terms happened to produce the same score.
+
+
+def _is_squared_function(term, function_type):
+    return (
+        term.is_Pow
+        and term.exp == 2
+        and isinstance(term.base, function_type)
+    )
+
+
+def _is_mixed_squared_base(base):
+    """Whether ``base`` is one trig square plus one hyperbolic square."""
+    if not base.is_Add or len(base.args) != 2:
+        return False
+    return (
+        any(
+            _is_squared_function(term, TrigonometricFunction)
+            for term in base.args
+        )
+        and any(
+            _is_squared_function(term, HyperbolicFunction)
+            for term in base.args
+        )
+    )
+
+
+def _has_expensive_fu_traversal(expr):
+    """Whether ``simplify`` is likely to hit SymPy's slow FU traversal."""
+    if _SYMPY_MAJOR_MINOR < (1, 13):
+        return False
+
+    return any(
+        (
+            node.is_Pow
+            and abs(node.exp) == sympy.S.Half
+            and _is_mixed_squared_base(node.base)
+        )
+        for node in preorder_traversal(expr)
+    )
+
+
+def simplify_for_display(expr):
+    """Simplify display output while avoiding a SymPy 1.13+ regression.
+
+    This helper is only for rendering. Algebraic operations retain ordinary
+    ``simplify``. Remove the fallback after SymPy replaces the nested
+    traversal introduced by gh-26390 and galgebra's minimum supported SymPy
+    includes that fix.
+    """
+    if _has_expensive_fu_traversal(expr):
+        return trigsimp(expr, method='old')
+    return simplify(expr)
diff --git a/galgebra/metric.py b/galgebra/metric.py
index aa26f91b..ecd5983e 100644
--- a/galgebra/metric.py
+++ b/galgebra/metric.py
@@ -14,6 +14,7 @@
 
 from . import printer
 from ._utils import cached_property as _cached_property
+from ._utils.simplify import simplify_for_display
 from .atoms import (
     BasisVectorSymbol, DotProductSymbol, MatrixFunction, Determinant,
 )
@@ -299,7 +300,9 @@ def symbols_list(s, indices=None, sub=True, commutative=False):
 
 
 class Simp:
-    modes = [simplify]
+    _default_modes = (simplify,)
+    modes = list(_default_modes)
+    _default_modes_instance = modes
 
     @staticmethod
     def profile(s):
@@ -312,6 +315,22 @@ def apply(expr):
             obj += apply_function_list(Simp.modes, coef) * base
         return obj
 
+    @staticmethod
+    def apply_display(expr):
+        """Apply the display fallback unless the user selected a profile."""
+        modes = (
+            [simplify_for_display]
+            if (
+                Simp.modes is Simp._default_modes_instance
+                and tuple(Simp.modes) == Simp._default_modes
+            )
+            else Simp.modes
+        )
+        obj = S.Zero
+        for coef, base in linear_expand_terms(expr):
+            obj += apply_function_list(modes, coef) * base
+        return obj
+
     @staticmethod
     def applymv(mv):
         return Mv(Simp.apply(mv.obj), ga=mv.Ga)
diff --git a/galgebra/mv.py b/galgebra/mv.py
index 544d8d5d..f1501802 100644
--- a/galgebra/mv.py
+++ b/galgebra/mv.py
@@ -612,7 +612,7 @@ def _sympystr(self, print_obj: printer.GaPrinter) -> str:
 
         # note: this just replaces `self` for the rest of this function
         obj = expand(self.obj)
-        obj = metric.Simp.apply(obj)
+        obj = metric.Simp.apply_display(obj)
         self = Mv(obj, ga=self.Ga)
 
         if self.i_grade == 0:
@@ -697,7 +697,7 @@ def append_plus(c_str):
         # note: this just replaces `self` for the rest of this function
         obj = expand(self.obj)
         try:
-            obj = metric.Simp.apply(obj)
+            obj = metric.Simp.apply_display(obj)
         except ZeroDivisionError:
             pass  # SymPy trigsimp regression; display without simplification
         self = Mv(obj, ga=self.Ga)
diff --git a/test/test_simplify.py b/test/test_simplify.py
new file mode 100644
index 00000000..abc3f11f
--- /dev/null
+++ b/test/test_simplify.py
@@ -0,0 +1,254 @@
+from unittest import mock
+
+import pytest
+from sympy import Add, Rational, cos, cosh, simplify, sin, sinh, sqrt, symbols
+
+from galgebra._utils import simplify as simplify_module
+from galgebra.ga import Ga
+from galgebra.metric import Simp
+
+
+x, y, u, v = symbols('x y u v')
+z = symbols('z:4')
+
+
+def _paired_trig_expression(count, extra=0):
+    terms = [sin(x + i) + cos(x + i) for i in range(count)]
+    if extra != 0:
+        terms.append(extra)
+    return Add(*terms, evaluate=False)
+
+
+def _mixed_nested_expression():
+    return 1/sqrt(sin(x)**2 + sinh(y)**2)
+
+
+def test_major_minor():
+    assert simplify_module._major_minor('1.13.3') == (1, 13)
+    assert simplify_module._major_minor('1.15.dev') == (1, 15)
+    assert simplify_module._major_minor('unknown') == (0, 0)
+
+
+def test_routes_only_observed_mixed_squared_shape():
+    plain = sqrt(sin(x) + sinh(y))
+    observed = _mixed_nested_expression()
+    other_power = (sin(x)**2 + sinh(y)**2)**Rational(1, 3)
+
+    with mock.patch.object(
+        simplify_module, '_SYMPY_MAJOR_MINOR', (1, 13)
+    ):
+        assert not simplify_module._has_expensive_fu_traversal(plain)
+        assert not simplify_module._has_expensive_fu_traversal(other_power)
+        assert simplify_module._has_expensive_fu_traversal(observed)
+
+
+def test_shallow_trig_sum_uses_real_general_simplifier():
+    rational = (y**2 - 1)/(y - 1)
+    expr = _paired_trig_expression(16, rational)
+
+    result = simplify_module.simplify_for_display(expr)
+
+    assert not result.has(rational)
+    assert simplify(result - expr) == 0
+
+
+def test_small_mixed_radical_cannot_borrow_unrelated_expression_cost():
+    rational = (y**2 - 1)/(y - 1)
+    expr = Add(
+        *[sin(x + i) + cos(x + i) for i in range(16)],
+        rational,
+        sqrt(sin(u) + sinh(v)),
+        evaluate=False,
+    )
+
+    with mock.patch.object(
+        simplify_module, '_SYMPY_MAJOR_MINOR', (1, 13)
+    ):
+        assert not simplify_module._has_expensive_fu_traversal(expr)
+
+    result = simplify_module.simplify_for_display(expr)
+
+    assert not result.has(rational)
+    assert simplify(result - expr) == 0
+
+
+def test_large_benign_mixed_radical_uses_general_simplifier():
+    rational = (y**2 - 1)/(y - 1)
+    expr = Add(
+        *[sin(x + i) + cos(x + i) for i in range(16)],
+        rational,
+        sqrt(sin(u) + sinh(v) + sum(z)),
+        evaluate=False,
+    )
+
+    with mock.patch.object(
+        simplify_module, '_SYMPY_MAJOR_MINOR', (1, 13)
+    ):
+        assert not simplify_module._has_expensive_fu_traversal(expr)
+
+    result = simplify_module.simplify_for_display(expr)
+
+    assert not result.has(rational)
+    assert simplify(result - expr) == 0
+
+
+def test_shallow_mixed_sum_does_not_match_failure_shape():
+    terms = [
+        sin(x + i) + cos(x + i) + sinh(y + i) + cosh(y + i)
+        for i in range(8)
+    ]
+    expr = Add(*terms, evaluate=False)
+
+    with mock.patch.object(
+        simplify_module, '_SYMPY_MAJOR_MINOR', (1, 13)
+    ):
+        assert not simplify_module._has_expensive_fu_traversal(expr)
+
+
+def test_algebra_keeps_general_simplification():
+    rational = (y**2 - 1)/(y - 1)
+
+    with mock.patch(
+        'galgebra.metric.simplify_for_display'
+    ) as display:
+        result = Simp.apply(rational)
+
+    assert not result.has(rational)
+    assert simplify(result - rational) == 0
+    display.assert_not_called()
+
+
+def test_display_route_can_preserve_unrelated_algebraic_form():
+    rational = (y**2 - 1)/(y - 1)
+    expr = _mixed_nested_expression() + rational
+
+    with (
+        mock.patch.object(
+            simplify_module, '_SYMPY_MAJOR_MINOR', (1, 13)
+        ),
+        mock.patch.object(simplify_module, 'simplify') as general,
+        mock.patch.object(
+            simplify_module, 'trigsimp', return_value=expr
+        ) as old,
+    ):
+        assert Simp.apply_display(expr) == expr
+
+    general.assert_not_called()
+    old.assert_called_once()
+    routed = old.call_args.args[0]
+    assert old.call_args.kwargs == {'method': 'old'}
+    assert simplify(routed - expr) == 0
+
+
+def test_small_display_expression_uses_general_simplification():
+    expr = sin(x)**2 + cos(x)**2
+
+    with (
+        mock.patch.object(
+            simplify_module, 'simplify', return_value=1
+        ) as general,
+        mock.patch.object(simplify_module, 'trigsimp') as old,
+    ):
+        assert Simp.apply_display(expr) == 1
+
+    general.assert_called_once()
+    assert simplify(general.call_args.args[0] - expr) == 0
+    old.assert_not_called()
+
+
+def test_sympy_before_1_13_uses_general_display_simplification():
+    expr = _paired_trig_expression(16)
+
+    with (
+        mock.patch.object(simplify_module, '_SYMPY_MAJOR_MINOR', (1, 12)),
+        mock.patch.object(
+            simplify_module, 'simplify', return_value=1
+        ) as general,
+        mock.patch.object(simplify_module, 'trigsimp') as old,
+    ):
+        assert Simp.apply_display(expr) == 1
+
+    general.assert_called_once()
+    assert simplify(general.call_args.args[0] - expr) == 0
+    old.assert_not_called()
+
+
+def test_custom_profile_overrides_display_fallback():
+    original_modes = Simp.modes
+    custom = mock.Mock(return_value=x)
+    Simp.profile([custom])
+    try:
+        assert Simp.apply_display(_paired_trig_expression(16)) == x
+    finally:
+        Simp.modes = original_modes
+
+    custom.assert_called_once()
+
+
+def test_default_profile_object_restores_display_fallback():
+    original_modes = Simp.modes
+    try:
+        Simp.profile([mock.Mock(return_value=x)])
+        Simp.profile(original_modes)
+        with mock.patch(
+            'galgebra.metric.simplify_for_display', return_value=1
+        ) as display:
+            assert Simp.apply_display(x) == 1
+    finally:
+        Simp.modes = original_modes
+
+    display.assert_called_once_with(x)
+
+
+def test_explicit_simplify_profile_overrides_display_fallback():
+    original_modes = Simp.modes
+    try:
+        Simp.profile([simplify])
+        with mock.patch(
+            'galgebra.metric.simplify_for_display'
+        ) as display:
+            assert Simp.apply_display(sin(x)**2 + cos(x)**2) == 1
+    finally:
+        Simp.modes = original_modes
+
+    display.assert_not_called()
+
+
+def test_in_place_profile_change_overrides_display_fallback():
+    custom = mock.Mock(return_value=x)
+    Simp.modes.append(custom)
+    try:
+        with mock.patch(
+            'galgebra.metric.simplify_for_display'
+        ) as display:
+            assert Simp.apply_display(x) == x
+    finally:
+        Simp.modes.remove(custom)
+
+    custom.assert_called_once_with(x)
+    display.assert_not_called()
+
+
+def test_prolate_spheroidal_divergence_renders():
+    if simplify_module._SYMPY_MAJOR_MINOR < (1, 13):
+        pytest.skip('display fallback targets SymPy 1.13 and newer')
+
+    a = symbols('a', real=True)
+    coords = xi, eta, phi = symbols('xi eta phi', real=True)
+    ps3d, *_ = Ga.build(
+        'e_xi e_eta e_phi',
+        X=[
+            a*sinh(xi)*sin(eta)*cos(phi),
+            a*sinh(xi)*sin(eta)*sin(phi),
+            a*cosh(xi)*cos(eta),
+        ],
+        coords=coords,
+        norm=True,
+    )
+    vector = ps3d.mv('A', 'vector', f=True)
+
+    rendered = str(ps3d.grad | vector)
+
+    assert 'D{eta}A__eta' in rendered
+    assert 'D{phi}A__phi' in rendered
+    assert 'D{xi}A__xi' in rendered