Repository: OpenMDAO/dymos

Issue:
Move timeseries expressions to ODE
### Proposed feature.

TImeseries expressions are useful but a more streamlined way to implement them is to include the ExecComp in a group with the ODE.

Once the Radau and GaussLobatto refactors are complete this should be possible and allow for the removal of a good bit of the complexity in timeseries.

### Example

There will be no change in the user-facing API.

Discussion:

Solution patch:
diff --git a/docs/dymos_book/_toc.yml b/docs/dymos_book/_toc.yml
index 8a39840341..78f63eb771 100644
--- a/docs/dymos_book/_toc.yml
+++ b/docs/dymos_book/_toc.yml
@@ -54,6 +54,7 @@ parts:
     - file: features/phases/constraints
     - file: features/phases/objective
     - file: features/phases/timeseries
+    - file: features/phases/calc_expressions
   - file: features/trajectories/trajectories
   - file: features/exploiting_sparsity
   - file: features/plotting
diff --git a/docs/dymos_book/api/phase_api.ipynb b/docs/dymos_book/api/phase_api.ipynb
index 86562c1a29..def47bb4c4 100644
--- a/docs/dymos_book/api/phase_api.ipynb
+++ b/docs/dymos_book/api/phase_api.ipynb
@@ -256,6 +256,12 @@
     "        :noindex:\n",
     "```\n",
     "\n",
+    "## add_calc_expr\n",
+    "```{eval-rst}\n",
+    "    .. automethod:: dymos.Phase.add_calc_expr\n",
+    "        :noindex:\n",
+    "```\n",
+    "\n",
     "## add_timeseries_output\n",
     "```{eval-rst}\n",
     "    .. automethod:: dymos.Phase.add_timeseries_output\n",
@@ -298,13 +304,6 @@
     "        :noindex:\n",
     "```"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/docs/dymos_book/examples/hull/hull_problem.ipynb b/docs/dymos_book/examples/hull/hull_problem.ipynb
index 2dd17f4532..e3a80a901f 100644
--- a/docs/dymos_book/examples/hull/hull_problem.ipynb
+++ b/docs/dymos_book/examples/hull/hull_problem.ipynb
@@ -508,7 +508,7 @@
    "provenance": []
   },
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "py311_2",
    "language": "python",
    "name": "python3"
   },
@@ -522,7 +522,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,
diff --git a/docs/dymos_book/features/phases/calc_expressions.ipynb b/docs/dymos_book/features/phases/calc_expressions.ipynb
new file mode 100644
index 0000000000..ffc9b55e14
--- /dev/null
+++ b/docs/dymos_book/features/phases/calc_expressions.ipynb
@@ -0,0 +1,191 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# This cell is mandatory in all Dymos documentation notebooks.\n",
+    "missing_packages = []\n",
+    "try:\n",
+    "    import openmdao.api as om  # noqa: F401\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !python -m pip install openmdao[notebooks]\n",
+    "    else:\n",
+    "        missing_packages.append('openmdao')\n",
+    "try:\n",
+    "    import dymos as dm  # noqa: F401\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !python -m pip install dymos\n",
+    "    else:\n",
+    "        missing_packages.append('dymos')\n",
+    "try:\n",
+    "    import pyoptsparse  # noqa: F401\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !pip install -q condacolab\n",
+    "        import condacolab\n",
+    "        condacolab.install_miniconda()\n",
+    "        !conda install -c conda-forge pyoptsparse\n",
+    "    else:\n",
+    "        missing_packages.append('pyoptsparse')\n",
+    "if missing_packages:\n",
+    "    raise EnvironmentError('This notebook requires the following packages '\n",
+    "                           'please install them and restart this notebook\\'s runtime: {\",\".join(missing_packages)}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Calculation Expressions\n",
+    "\n",
+    "Sometimes a user may wish to calculate some additional output that is not inherently provided by the ODE system. In cases where they have access to the model, it might be trivial to add it to the ODE's outputs. In other cases, they might not want to modify a model provided by someone else.\n",
+    "\n",
+    "Calculation expressions provide the user with a way of computing additional outputs that is extremely general.\n",
+    "Older version of dymos supported expressions for use as constraints, objectives, or timeseries outputs.\n",
+    "Now dymos just computes these expressions along-side the ODE, making them available to basically all of dymos (constraints, objectives, timeseries outputs, state rates, etc.)\n",
+    "\n",
+    "The most general way to add a calculated expression to an ODE output is to use the `phase.add_calc_expr` method.\n",
+    "\n",
+    "```{eval-rst}\n",
+    "    .. automethod:: dymos.Phase.add_calc_expr\n",
+    "        :noindex:\n",
+    "```\n",
+    "\n",
+    "Providing the calculated expression as the name in `add_boundary_constraint`, `add_path_constraint`, `add_objective`, or `add_timeseries_output` will still work for backward compatibility.\n",
+    "In general, using `add_calc_expr` will just always work, and makes that output variable available across the phase.\n",
+    "\n",
+    "## Units\n",
+    "\n",
+    "In general, a given calculation that a user might apply assumes specific units for each of the input variables.\n",
+    "In dymos, introspection for units happens _after_ the ODE is completely defined.\n",
+    "Because the ExecComp that computes the calcuation expressions is itself part of the ODE, we cannot wait for this introspection step.\n",
+    "Units for inputs **will default to None** if not otherwise specified.\n",
+    "\n",
+    "If you're not certain of the units of an input variable upon which your calculation relies, you should explicitly provide it via keyword arguments to the `add_calc_expr` call (or one of the associated calls that accepts calculation expressions, mentioned above).\n",
+    "\n",
+    "## Example - Adding an additional calculation to the brachistochrone\n",
+    "\n",
+    "Let's assume that we have a brachistochrone ODE that was provided to us, but for configuration control reasons, we might not want to edit the source code.\n",
+    "\n",
+    "Here we're adding a calculation expression for the additional output k, which is the ratio\n",
+    "\n",
+    "\\begin{align}\n",
+    "    k &= \\frac{sin(\\{theta})}{v}\n",
+    "\\end{align}\n",
+    "\n",
+    "In the brachistochrone solution, $k$ is approximately constant, and its value is dependent upon the start and end point of the trajectory, as well as the gravitational parameter. In practice, it's value is somewhat large at the start since we're starting with nearly zero velocity."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import openmdao.api as om\n",
+    "import dymos as dm\n",
+    "from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE\n",
+    "\n",
+    "p = om.Problem(model=om.Group())\n",
+    "\n",
+    "p.driver = om.ScipyOptimizeDriver()\n",
+    "\n",
+    "p.driver.declare_coloring()\n",
+    "\n",
+    "# Create the transcription so we can get the number of nodes for the downstream analysis\n",
+    "tx = dm.Radau(num_segments=20, order=3, compressed=False)\n",
+    "\n",
+    "traj = dm.Trajectory()\n",
+    "phase = dm.Phase(transcription=tx, ode_class=BrachistochroneODE)\n",
+    "traj.add_phase('phase0', phase)\n",
+    "\n",
+    "p.model.add_subsystem('traj', traj)\n",
+    "\n",
+    "phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))\n",
+    "\n",
+    "phase.add_state('x', units='m', rate_source='xdot', fix_initial=True, fix_final=True)\n",
+    "phase.add_state('y', units='m', rate_source='ydot', fix_initial=True, fix_final=True)\n",
+    "phase.add_state('v', units='m/s', rate_source='vdot', fix_initial=True, fix_final=False)\n",
+    "\n",
+    "phase.add_control('theta', units='deg', lower=0.01, upper=179.9,\n",
+    "                  continuity=True, rate_continuity=True)\n",
+    "\n",
+    "phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)\n",
+    "\n",
+    "phase.add_calc_expr('k = sin(theta) / v', theta={'units': 'rad'}, v={'units': 'm/s'})\n",
+    "\n",
+    "# Minimize time at the end of the phase\n",
+    "phase.add_objective('time', loc='final', scaler=1)\n",
+    "\n",
+    "p.setup(check=True)\n",
+    "\n",
+    "phase.set_time_val(initial=0.0, duration=2.0)\n",
+    "phase.set_state_val('x', [0, 10])\n",
+    "phase.set_state_val('y', [10, 5])\n",
+    "phase.set_state_val('v', [1.0E-1, 9.9])\n",
+    "phase.set_control_val('theta', [5, 100])\n",
+    "phase.set_parameter_val('g', 9.80665)\n",
+    "\n",
+    "dm.run_problem(p, simulate=True, make_plots=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "hide-input"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from IPython.display import HTML\n",
+    "\n",
+    "# Define the path to the HTML file\n",
+    "html_file_path = p.get_reports_dir() / 'traj_results_report.html'\n",
+    "html_content = html_file_path.read_text()\n",
+    "\n",
+    "# Inject CSS to control the output cell height and avoid scrollbars\n",
+    "html_with_custom_height = f\"\"\"\n",
+    "<div style=\"height: 800px; overflow: auto;\">\n",
+    "    {html_content}\n",
+    "</div>\n",
+    "\"\"\"\n",
+    "\n",
+    "HTML(html_with_custom_height)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "py312",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.8"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/docs/dymos_book/features/phases/phases.ipynb b/docs/dymos_book/features/phases/phases.ipynb
index f0a1db8e4a..865d5a9318 100644
--- a/docs/dymos_book/features/phases/phases.ipynb
+++ b/docs/dymos_book/features/phases/phases.ipynb
@@ -114,7 +114,10 @@
     "\n",
     "[Objective](objective.ipynb)\n",
     "\n",
-    "[Timeseries Outputs](timeseries.ipynb)\n"
+    "[Timeseries Outputs](timeseries.ipynb)\n",
+    "\n",
+    "[Calculation Expressions](calc_expressions.ipynb)\n",
+    "\n"
    ]
   }
  ],
diff --git a/docs/dymos_book/features/phases/timeseries.ipynb b/docs/dymos_book/features/phases/timeseries.ipynb
index 14886f01b5..c6eed220da 100644
--- a/docs/dymos_book/features/phases/timeseries.ipynb
+++ b/docs/dymos_book/features/phases/timeseries.ipynb
@@ -269,7 +269,7 @@
    }
   },
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "py311_2",
    "language": "python",
    "name": "python3"
   },
diff --git a/dymos/phase/options.py b/dymos/phase/options.py
index 3a9ef71900..0a8510ff7c 100644
--- a/dymos/phase/options.py
+++ b/dymos/phase/options.py
@@ -633,10 +633,6 @@ def __init__(self, read_only=False):
                      desc='If True, the given indices will be treated as indices into a C-order flattened array based '
                           'on the shaped of the constrained variable at a point in time.')
 
-        self.declare(name='is_expr', types=bool, default=False,
-                     desc='If True, the given constraint is an expression that must be evaluated rather than a'
-                          ' single variable.')
-
 
 class TimeseriesOutputOptionsDictionary(om.OptionsDictionary):
     """
@@ -680,12 +676,6 @@ def __init__(self, read_only=False):
         self.declare(name='is_rate', default=False, allow_none=False,
                      desc='If True this is a rate.')
 
-        self.declare(name='is_expr', default=False, allow_none=False,
-                     desc='If true the requested timeseries is a mathematical expression')
-
-        self.declare(name='expr_kwargs', default={}, allow_none=False,
-                     desc='Options to be passed to the timeseries expression comp when adding the expression.')
-
 
 class PhaseTimeseriesOptionsDictionary(om.OptionsDictionary):
     """
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index e123f2b2f3..d7e28c50a5 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -26,7 +26,7 @@
 from ..utils.indexing import get_constraint_flat_idxs
 from ..utils.introspection import configure_time_introspection, _configure_constraint_introspection, \
     configure_controls_introspection, configure_parameters_introspection, \
-    configure_timeseries_output_introspection, classify_var, configure_timeseries_expr_introspection
+    configure_timeseries_output_introspection, classify_var
 from ..utils.misc import _unspecified, create_subprob, om_version
 from ..utils.lgl import lgl
 
@@ -85,6 +85,9 @@ def __init__(self, from_phase=None, **kwargs):
         self._objectives = {}
         self.sim_prob = None
 
+        # Expressions to be computed along with the ODE.
+        self._calc_exprs = {}
+
         super(Phase, self).__init__(**_kwargs)
 
     def duplicate(self, transcription=None, boundary_constraints=False, path_constraints=False, objectives=False,
@@ -1291,27 +1294,21 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
             raise ValueError(f'Invalid boundary constraint location "{loc}". Must be '
                              '"initial" or "final".')
 
-        expr_operators = ['(', '+', '-', '/', '*', '&', '%', '@']
         if '=' in name:
-            is_expr = True
-        elif '=' not in name and any(opr in name for opr in expr_operators):
-            raise ValueError(f'The expression provided `{name}` has invalid format. '
-                             'Expression may be a single variable or an equation '
-                             'of the form `constraint_name = func(vars)`')
-        else:
-            is_expr = False
-
-        if is_expr:
-            constraint_name = name.split('=')[0].strip()
+            _name = constraint_name = name.split('=')[0].strip()
+            self.add_calc_expr(name)
         elif constraint_name is None:
+            _name = name
             constraint_name = name.rpartition('.')[-1]
+        else:
+            _name = name
 
         bc_list = self._initial_boundary_constraints if loc == 'initial' else self._final_boundary_constraints
 
-        existing_bc = [bc for bc in bc_list if bc['name'] == name and bc['indices'] is None and indices is None]
+        existing_bc = [bc for bc in bc_list if bc['name'] == _name and bc['indices'] is None and indices is None]
 
         if existing_bc:
-            raise ValueError(f'Cannot add new {loc} boundary constraint for variable `{name}` and indices {indices}. '
+            raise ValueError(f'Cannot add new {loc} boundary constraint for variable `{_name}` and indices {indices}. '
                              f'One already exists.')
 
         existing_bc_name = [bc for bc in bc_list if bc['name'] == constraint_name and
@@ -1325,7 +1322,7 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
         bc = ConstraintOptionsDictionary()
         bc_list.append(bc)
 
-        bc['name'] = name
+        bc['name'] = _name
         bc['constraint_name'] = constraint_name
         bc['lower'] = lower
         bc['upper'] = upper
@@ -1339,7 +1336,6 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
         bc['linear'] = linear
         bc['units'] = units
         bc['flat_indices'] = flat_indices
-        bc['is_expr'] = is_expr
 
         # Automatically add the requested variable to the timeseries outputs if it's an ODE output.
         var_type = self.classify_var(name)
@@ -1404,26 +1400,20 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
             If True, treat indices as flattened C-ordered indices of elements to constrain at each given point in time.
             Otherwise, indices should be a tuple or list giving the elements to constrain at each point in time.
         """
-        expr_operators = ['(', '+', '-', '/', '*', '&', '%', '@']
         if '=' in name:
-            is_expr = True
-        elif '=' not in name and any(opr in name for opr in expr_operators):
-            raise ValueError(f'The expression provided `{name}` has invalid format. '
-                             'Expression may be a single variable or an equation '
-                             'of the form `constraint_name = func(vars)`')
-        else:
-            is_expr = False
-
-        if is_expr:
-            constraint_name = name.split('=')[0].strip()
+            _name = constraint_name = name.split('=')[0].strip()
+            self.add_calc_expr(name)
         elif constraint_name is None:
+            _name = name
             constraint_name = name.rpartition('.')[-1]
+        else:
+            _name = name
 
         existing_pc = [pc for pc in self._path_constraints
                        if pc['name'] == name and pc['indices'] == indices and pc['flat_indices'] == flat_indices]
 
         if existing_pc:
-            raise ValueError(f'Cannot add new path constraint for variable `{name}` and indices {indices}. '
+            raise ValueError(f'Cannot add new path constraint for variable `{_name}` and indices {indices}. '
                              f'One already exists.')
 
         existing_bc_name = [pc for pc in self._path_constraints
@@ -1437,7 +1427,7 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
         pc = ConstraintOptionsDictionary()
         self._path_constraints.append(pc)
 
-        pc['name'] = name
+        pc['name'] = _name
         pc['constraint_name'] = constraint_name
         pc['lower'] = lower
         pc['upper'] = upper
@@ -1451,7 +1441,6 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
         pc['linear'] = linear
         pc['units'] = units
         pc['flat_indices'] = flat_indices
-        pc['is_expr'] = is_expr
 
         # Automatically add the requested variable to the timeseries outputs if it's an ODE output.
         var_type = self.classify_var(name)
@@ -1496,7 +1485,6 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap
         """
         if type(name) is list:
             for i, name_i in enumerate(name):
-                expr = True if '=' in name_i else False
                 if type(units) is dict:  # accept dict for units when using array of name
                     unit = units.get(name_i, None)
                 elif type(units) is list:  # allow matching list for units
@@ -1504,26 +1492,40 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap
                 else:
                     unit = units
 
-                oname = self._add_timeseries_output(name_i, output_name=output_name,
-                                                    units=unit,
-                                                    shape=shape,
-                                                    timeseries=timeseries,
-                                                    rate=False,
-                                                    expr=expr)
+                if '=' in name_i:
+                    output_i = name.split('=')[0].strip()
+                    self.add_calc_expr(name, add_timeseries=True, units=unit, shape=shape, **kwargs)
+                    oname = self._add_timeseries_output(output_i,
+                                                        units=unit,
+                                                        shape=shape,
+                                                        timeseries=timeseries,
+                                                        rate=False)
+                else:
+                    oname = self._add_timeseries_output(name_i, output_name=output_name,
+                                                        units=unit,
+                                                        shape=shape,
+                                                        timeseries=timeseries,
+                                                        rate=False)
 
                 # Handle specific units for wildcard names.
                 if oname is not None and '*' in name_i:
                     self._timeseries[timeseries]['outputs'][oname]['wildcard_units'] = units
 
         else:
-            expr = True if '=' in name else False
-            self._add_timeseries_output(name, output_name=output_name,
+            if '=' in name:
+                output = name.split('=')[0].strip()
+                _kwargs = {k: v for k, v in kwargs.items()}
+                if units is not _unspecified:
+                    if output not in _kwargs:
+                        _kwargs[output] = {'units': units}
+                self.add_calc_expr(name, add_timeseries=False, **_kwargs)
+            else:
+                output = name
+            self._add_timeseries_output(output, output_name=output_name,
                                         units=units,
                                         shape=shape,
                                         timeseries=timeseries,
-                                        rate=False,
-                                        expr=expr,
-                                        expr_kwargs=kwargs)
+                                        rate=False)
 
     def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified, shape=_unspecified,
                                    timeseries='timeseries'):
@@ -1581,7 +1583,7 @@ def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified,
                                         rate=True)
 
     def _add_timeseries_output(self, name, output_name=None, units=_unspecified, shape=_unspecified,
-                               timeseries='timeseries', rate=False, expr=False, expr_kwargs=None):
+                               timeseries='timeseries', rate=False):
         r"""
         Add a single variable or rate to the timeseries outputs of the phase.
 
@@ -1611,11 +1613,6 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha
         rate : bool
             If True, add the rate of change of the named variable to the timeseries outputs of the
             phase.  The rate variable will be named f'{name}_rate'.  Defaults to False.
-        expr : bool
-            True if the given name is an expression for an ExecComp.
-        expr_kwargs : dict
-            Keyword arguments for the expression.
-
 
         Returns
         -------
@@ -1626,9 +1623,7 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha
             raise ValueError(f'Timeseries {timeseries} does not exist in phase {self.pathname}')
 
         if output_name is None:
-            if expr:
-                output_name = name.split('=')[0].strip()
-            elif '*' in name:
+            if '*' in name:
                 output_name = name
             elif output_name is None:
                 output_name = name.rpartition('.')[-1]
@@ -1644,8 +1639,6 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha
             ts_output['units'] = units
             ts_output['shape'] = shape
             ts_output['is_rate'] = rate
-            ts_output['is_expr'] = expr
-            ts_output['expr_kwargs'] = expr_kwargs
 
             self._timeseries[timeseries]['outputs'][output_name] = ts_output
 
@@ -1709,17 +1702,11 @@ def add_objective(self, name, loc='final', index=None, shape=(1,), units=None, r
             If specified, this design var will be grouped for parallel derivative
             calculations with other variables sharing the same parallel_deriv_color.
         """
-        expr_operators = ['(', '+', '-', '/', '*', '&', '%', '@']
         if '=' in name:
-            is_expr = True
-        elif '=' not in name and any(opr in name for opr in expr_operators):
-            raise ValueError(f'The expression provided `{name}` has invalid format. '
-                             'Expression may be a single variable or an equation '
-                             'of the form `constraint_name = func(vars)`')
+            obj_name = name.split('=')[0].strip()
+            self.add_calc_expr(name)
         else:
-            is_expr = False
-
-        obj_name = name.split('=')[0].strip() if is_expr else name
+            obj_name = name
 
         obj_dict = {'name': name,
                     'loc': loc,
@@ -1730,11 +1717,47 @@ def add_objective(self, name, loc='final', index=None, shape=(1,), units=None, r
                     'ref0': ref0,
                     'adder': adder,
                     'scaler': scaler,
-                    'parallel_deriv_color': parallel_deriv_color,
-                    'is_expr': is_expr}
+                    'parallel_deriv_color': parallel_deriv_color}
         self._objectives[obj_name] = obj_dict
-        if is_expr and obj_name not in self._timeseries['timeseries']['outputs']:
-            self.add_timeseries_output(name, output_name=obj_name, units=units, shape=shape)
+
+    def add_calc_expr(self, expr, add_timeseries=True, **kwargs):
+        """
+        Adds an expression to be computed immediately after the user-given ODE.
+
+        Internally, dymos will wrap the user-given ODE in a group along with
+        an ExecComp that evalutes the expressions given.
+
+        Unlike standard OpenMDAO ExecComp expressions, those specified here
+        may use the ODE-relative path of variables (if they are not promoted
+        to the top of the ODE).
+
+        This function may be called more than one time with the same expression. In that case, any
+        new info from kwargs will be updated.
+
+        Parameters
+        ----------
+        expr : str
+            The expression to be computed.
+        add_timeseries : bool
+            If True, add the output of the expression to the timeseries.
+        **kwargs : dict
+            Any arguments to be forwarded to the ExecComp when the expression in added.
+        """
+        if expr in self._calc_exprs:
+            self._calc_exprs[expr].update(kwargs)
+        else:
+            self._calc_exprs[expr] = kwargs
+
+        output_name = expr.split('=')[0].strip()
+        if add_timeseries and output_name not in self._timeseries['timeseries']['outputs']:
+            output = expr.split('=')[0].strip()
+            if 'units' in kwargs:
+                units = kwargs['units']
+            elif output in kwargs:
+                units = kwargs[output]['units']
+            else:
+                units = None
+            self.add_timeseries_output(output_name, units=units)
 
     def set_time_options(self, units=_unspecified, fix_initial=_unspecified,
                          fix_duration=_unspecified, input_initial=_unspecified,
@@ -2103,19 +2126,12 @@ def set_duration_balance(self, name, val=0.0, index=None, units=None, mult_val=N
                    'mult_val': mult_val,
                    'normalize': normalize}
 
-        expr_operators = ['(', '+', '-', '/', '*', '&', '%', '@']
         if '=' in name:
-            is_expr = True
-        elif '=' not in name and any(opr in name for opr in expr_operators):
-            raise ValueError(f'The expression provided `{name}` has invalid format. '
-                             'Expression may be a single variable or an equation '
-                             'of the form `constraint_name = func(vars)`')
+            balance_name = name.split('=')[0].strip()
+            self.add_calc_expr(name)
         else:
-            is_expr = False
-
-        balance_name = name.split('=')[0].strip() if is_expr else name
+            balance_name = name
 
-        options['is_expr'] = is_expr
         options['balance_name'] = balance_name
 
         self.time_options['t_duration_balance_options'] = options
@@ -2239,8 +2255,6 @@ def configure(self):
 
         _configure_constraint_introspection(self)
 
-        configure_timeseries_expr_introspection(self)
-
         transcription.configure_boundary_constraints(self)
 
         transcription.configure_path_constraints(self)
@@ -2581,6 +2595,8 @@ def get_simulation_phase(self, times_per_seg=_unspecified, method=_unspecified,
         sim_phase.simulate_options = deepcopy(self.simulate_options)
         sim_phase.timeseries_options = deepcopy(self.timeseries_options)
 
+        sim_phase._calc_exprs = deepcopy(self._calc_exprs)
+
         return sim_phase
 
     def initialize_values_from_phase(self, prob, from_phase, phase_path='', skip_params=None):
@@ -2621,13 +2637,13 @@ def initialize_values_from_phase(self, prob, from_phase, phase_path='', skip_par
 
         # Set the integration times
         time_name = phs.time_options['name']
-        op = op_dict[f'timeseries.timeseries_comp.{time_name}']
+        op = op_dict[f'timeseries.{time_name}']
         prob.set_val(f'{self_path}t_initial', op['val'][0, ...])
         prob.set_val(f'{self_path}t_duration', op['val'][-1, ...] - op['val'][0, ...])
 
         # Assign initial state values
         for name in phs.state_options:
-            op = op_dict[f'timeseries.timeseries_comp.states:{name}']
+            op = op_dict[f'timeseries.states:{name}']
             prob[f'{self_path}initial_states:{name}'][...] = op['val'][0, ...]
 
         # Assign control values
@@ -3007,7 +3023,7 @@ def load_case(self, case):
         integration_name = self.time_options['name']
 
         try:
-            prev_time_path = prev_vars_abs2prom[f'{self.pathname}.timeseries.timeseries_comp.{integration_name}']
+            prev_time_path = prev_vars_abs2prom[f'{self.pathname}.timeseries.{integration_name}']
         except KeyError:
             om.issue_warning(f'load_case for phase {self.name} failed - phase not found in case data.')
             return
diff --git a/dymos/transcriptions/analytic/analytic.py b/dymos/transcriptions/analytic/analytic.py
index d063d05d78..0ee2b43397 100644
--- a/dymos/transcriptions/analytic/analytic.py
+++ b/dymos/transcriptions/analytic/analytic.py
@@ -5,8 +5,9 @@
 from ...utils.introspection import configure_analytic_states_introspection, get_promoted_vars, \
     get_source_metadata, configure_analytic_states_discovery
 from ...utils.indexing import get_src_indices_by_row
+from ...utils.ode_utils import _make_ode_system
 from ..grid_data import GridData
-from ..common import TimeComp, TimeseriesOutputGroup, TimeseriesOutputComp
+from ..common import TimeComp, TimeseriesOutputComp
 
 
 class Analytic(TranscriptionBase):
@@ -206,13 +207,15 @@ def setup_ode(self, phase):
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
         """
-        ODEClass = phase.options['ode_class']
         grid_data = self.grid_data
 
-        kwargs = phase.options['ode_init_kwargs']
-        phase.add_subsystem('rhs',
-                            subsys=ODEClass(num_nodes=grid_data.subset_num_nodes['all'],
-                                            **kwargs))
+        rhs = _make_ode_system(ode_class=phase.options['ode_class'],
+                               num_nodes=grid_data.subset_num_nodes['all'],
+                               ode_init_kwargs=phase.options['ode_init_kwargs'],
+                               calc_exprs=phase._calc_exprs,
+                               parameter_options=phase.parameter_options)
+
+        phase.add_subsystem('rhs', rhs)
 
     def configure_ode(self, phase):
         """
@@ -292,11 +295,6 @@ def setup_timeseries_outputs(self, phase):
         gd = self.grid_data
 
         for name, options in phase._timeseries.items():
-            has_expr = False
-            for _, output_options in options['outputs'].items():
-                if output_options['is_expr']:
-                    has_expr = True
-                    break
             if options['transcription'] is None:
                 ogd = None
             else:
@@ -307,8 +305,7 @@ def setup_timeseries_outputs(self, phase):
                                                    output_subset=options['subset'],
                                                    time_units=phase.time_options['units'])
 
-            timeseries_group = TimeseriesOutputGroup(has_expr=has_expr, timeseries_output_comp=timeseries_comp)
-            phase.add_subsystem(name, subsys=timeseries_group)
+            phase.add_subsystem(name, subsys=timeseries_comp)
 
             phase.connect('dt_dstau', f'{name}.dt_dstau', flat_src_indices=True)
 
diff --git a/dymos/transcriptions/common/__init__.py b/dymos/transcriptions/common/__init__.py
index da70e7436f..cbece8a0e2 100644
--- a/dymos/transcriptions/common/__init__.py
+++ b/dymos/transcriptions/common/__init__.py
@@ -2,5 +2,4 @@
 from .control_comp import ControlInterpComp
 from .parameter_comp import ParameterComp
 from .time_comp import TimeComp
-from .timeseries_group import TimeseriesOutputGroup
 from .timeseries_output_comp import TimeseriesOutputComp
diff --git a/dymos/transcriptions/common/timeseries_group.py b/dymos/transcriptions/common/timeseries_group.py
deleted file mode 100644
index d124301308..0000000000
--- a/dymos/transcriptions/common/timeseries_group.py
+++ /dev/null
@@ -1,38 +0,0 @@
-import openmdao.api as om
-
-from .timeseries_output_comp import TimeseriesOutputComp
-
-
-class TimeseriesOutputGroup(om.Group):
-    """
-    Class definition of the TimeComp.
-
-    Parameters
-    ----------
-    **kwargs : dict
-        Dictionary of optional arguments.
-    """
-    def __init__(self, **kwargs):
-        super().__init__(**kwargs)
-
-    def initialize(self):
-        """
-        Declare component options.
-        """
-        self.options.declare('timeseries_output_comp', types=TimeseriesOutputComp, recordable=False,
-                             desc='Timeseries component specific to the transcription of the optimal control problem.')
-
-        self.options.declare('has_expr', types=bool,
-                             desc='If true, timeseries group has an expression to be computed')
-
-    def setup(self):
-        """
-        Define the structure of the timeseries group.
-        """
-        timeseries_output_comp = self.options['timeseries_output_comp']
-        has_expr = self.options['has_expr']
-        if has_expr:
-            self.add_subsystem('timeseries_exec_comp', om.ExecComp(has_diag_partials=True))
-
-        self.add_subsystem('timeseries_comp', timeseries_output_comp,
-                           promotes_inputs=['*'], promotes_outputs=['*'])
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index ee834ce034..c6de23f010 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -5,7 +5,6 @@
 import openmdao.api as om
 from openmdao.utils.om_warnings import warn_deprecation
 
-from ..common.timeseries_output_comp import TimeseriesOutputComp
 from .explicit_shooting_continuity_comp import ExplicitShootingContinuityComp
 from ..transcription_base import TranscriptionBase
 from ..grid_data import BirkhoffGrid, GaussLobattoGrid, RadauGrid, UniformGrid
@@ -14,7 +13,9 @@
 from ...utils.indexing import get_src_indices_by_row
 from ...utils.introspection import get_promoted_vars, get_source_metadata, get_targets, _get_targets_metadata
 from ...utils.constants import INF_BOUND
-from ..common import TimeComp, TimeseriesOutputGroup, ControlInterpComp, ParameterComp
+
+from ...utils.ode_utils import _make_ode_system
+from ..common import TimeComp, TimeseriesOutputComp, ControlInterpComp, ParameterComp
 
 
 class ExplicitShooting(TranscriptionBase):
@@ -308,10 +309,17 @@ def setup_ode(self, phase):
                                    ode_init_kwargs=phase.options['ode_init_kwargs'],
                                    standalone_mode=False,
                                    reports=self.options['subprob_reports'],
-                                   control_interp=self.options['control_interp'])
+                                   control_interp=self.options['control_interp'],
+                                   calc_exprs=phase._calc_exprs)
         phase.add_subsystem('integrator', integ)
-        phase.add_subsystem('ode', phase.options['ode_class'](num_nodes=self._output_grid_data.num_nodes,
-                                                              **phase.options['ode_init_kwargs']))
+
+        ode = _make_ode_system(ode_class=phase.options['ode_class'],
+                               num_nodes=self._output_grid_data.num_nodes,
+                               calc_exprs=phase._calc_exprs,
+                               ode_init_kwargs=phase.options['ode_init_kwargs'],
+                               parameter_options=phase.parameter_options)
+
+        phase.add_subsystem('ode', ode)
 
     def configure_ode(self, phase):
         """
@@ -523,18 +531,11 @@ def setup_timeseries_outputs(self, phase):
             The phase object to which this transcription instance applies.
         """
         for name, options in phase._timeseries.items():
-            has_expr = False
-            for _, output_options in options['outputs'].items():
-                if output_options['is_expr']:
-                    has_expr = True
-                    break
-
             timeseries_comp = TimeseriesOutputComp(input_grid_data=self._output_grid_data,
                                                    output_grid_data=self._output_grid_data,
                                                    output_subset=options['subset'],
                                                    time_units=phase.time_options['units'])
-            timeseries_group = TimeseriesOutputGroup(has_expr=has_expr, timeseries_output_comp=timeseries_comp)
-            phase.add_subsystem(name, subsys=timeseries_group)
+            phase.add_subsystem(name, subsys=timeseries_comp)
 
             phase.connect('dt_dstau', f'{name}.dt_dstau', flat_src_indices=True)
 
diff --git a/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py b/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py
index 2aa64ad0ee..acf93d304d 100644
--- a/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py
+++ b/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py
@@ -14,6 +14,7 @@
     configure_states_discovery, configure_states_introspection, _get_targets_metadata, \
     _get_common_metadata, get_promoted_vars
 from ...utils.misc import get_rate_units, _unspecified, _none_or_unspecified
+from ...utils.ode_utils import _make_ode_system
 
 
 class ODEEvaluationGroup(om.Group):
@@ -44,13 +45,15 @@ class ODEEvaluationGroup(om.Group):
         The number of points at which the ODE is simultaneously evaluated.
     control_interp : str
         The control interpolation technique to be used. Must be either 'vandermonde' or 'barycentric'.
+    calc_exprs : dict
+        A dictionary of ODE expressions.
     **kwargs : dict
         Additional keyword arguments passed to Group.
     """
 
     def __init__(self, ode_class, input_grid_data, time_options, state_options, parameter_options, control_options,
                  ode_init_kwargs=None, compute_derivs=True, vec_size=1,
-                 control_interp='vandermonde', **kwargs):
+                 control_interp='vandermonde', calc_exprs=None, **kwargs):
         super().__init__(**kwargs)
 
         # This component creates copies of the variable options from the phase.
@@ -68,6 +71,7 @@ def __init__(self, ode_class, input_grid_data, time_options, state_options, para
         self._compute_derivs = compute_derivs
         self._vec_size = vec_size
         self._ode_init_kwargs = {} if ode_init_kwargs is None else ode_init_kwargs
+        self._calc_exprs = {} if calc_exprs is None else calc_exprs
         self._control_interp = control_interp
 
     def set_segment_index(self, seg_idx):
@@ -130,7 +134,13 @@ def setup(self):
                                                                                      compute_derivs=self._compute_derivs),
                                                         promotes_inputs=['ptau', 'stau', 'dstau_dt', 't_duration'])
 
-        self.add_subsystem('ode', self._ode_class(num_nodes=self._vec_size, **self._ode_init_kwargs))
+        ode = _make_ode_system(ode_class=self._ode_class,
+                               num_nodes=self._vec_size,
+                               ode_init_kwargs=self._ode_init_kwargs,
+                               calc_exprs=self._calc_exprs,
+                               parameter_options=self._parameter_options)
+
+        self.add_subsystem('ode', ode)
 
         self.add_subsystem('state_rate_collector',
                            StateRateCollectorComp(state_options=self._state_options,
diff --git a/dymos/transcriptions/explicit_shooting/ode_integration_comp.py b/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
index bb11a748cb..bf52150cc5 100644
--- a/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
+++ b/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
@@ -38,6 +38,8 @@ class ODEIntegrationComp(om.ExplicitComponent):
         If True, assume this component is being run as its own system. As part of the Dymos
         ShootingPhase, this system needs to be setup during the Phase configure process and setting this
         to False will enable that behavior.
+    calc_exprs : dict
+        ODE Expressions provided by the controlling phase.
     **kwargs : dict
         Additional keyword arguments passed to Group.
 
@@ -47,7 +49,7 @@ class ODEIntegrationComp(om.ExplicitComponent):
     theta:  U+03B8
     """
     def __init__(self, input_grid_data, time_options, state_options, parameter_options=None, control_options=None,
-                 output_grid_data=None, reports=False, standalone_mode=True, **kwargs):
+                 output_grid_data=None, reports=False, standalone_mode=True, calc_exprs=None, **kwargs):
         super().__init__(**kwargs)
         self.time_options = time_options
         self.state_options = state_options
@@ -75,6 +77,8 @@ def __init__(self, input_grid_data, time_options, state_options, parameter_optio
         self._no_check_partials = not dymos_options['include_check_partials']
         self._num_control_input_nodes = input_grid_data.subset_num_nodes['control_input']
 
+        self._calc_exprs = {} if calc_exprs is None else calc_exprs
+
     def initialize(self):
         """
         Declare options for the ODEIntegrationComp.
@@ -108,7 +112,8 @@ def _setup_subprob(self):
                                                  ode_init_kwargs=self.options['ode_init_kwargs'],
                                                  input_grid_data=self._input_grid_data,
                                                  compute_derivs=self.options['propagate_derivs'],
-                                                 control_interp=self.options['control_interp']),
+                                                 control_interp=self.options['control_interp'],
+                                                 calc_exprs=self._calc_exprs),
                               promotes_inputs=['*'],
                               promotes_outputs=['*'])
 
diff --git a/dymos/transcriptions/pseudospectral/birkhoff.py b/dymos/transcriptions/pseudospectral/birkhoff.py
index e81cc03ea1..5617028c76 100644
--- a/dymos/transcriptions/pseudospectral/birkhoff.py
+++ b/dymos/transcriptions/pseudospectral/birkhoff.py
@@ -3,7 +3,7 @@
 import openmdao.api as om
 
 from ..transcription_base import TranscriptionBase
-from ..common import TimeComp, TimeseriesOutputGroup, TimeseriesOutputComp
+from ..common import TimeComp, TimeseriesOutputComp
 from .components import BirkhoffIterGroup, BirkhoffBoundaryGroup
 
 from ..grid_data import BirkhoffGrid
@@ -211,23 +211,22 @@ def setup_ode(self, phase):
         grid_data = self.grid_data
 
         ode_init_kwargs = phase.options['ode_init_kwargs']
-        ibcs = phase._initial_boundary_constraints
-        fbcs = phase._final_boundary_constraints
 
         phase.add_subsystem('ode_iter_group',
                             subsys=BirkhoffIterGroup(grid_data=grid_data, state_options=phase.state_options,
                                                      time_options=phase.time_options,
+                                                     parameter_options=phase.parameter_options,
                                                      ode_class=ODEClass,
-                                                     ode_init_kwargs=ode_init_kwargs),
+                                                     ode_init_kwargs=ode_init_kwargs,
+                                                     calc_exprs=phase._calc_exprs),
                             promotes=['*'])
 
         phase.add_subsystem('boundary_vals',
                             subsys=BirkhoffBoundaryGroup(grid_data=grid_data,
                                                          ode_class=ODEClass,
                                                          ode_init_kwargs=ode_init_kwargs,
-                                                         initial_boundary_constraints=ibcs,
-                                                         final_boundary_constraints=fbcs,
-                                                         objectives=phase._objectives),
+                                                         calc_exprs=phase._calc_exprs,
+                                                         parameter_options=phase.parameter_options),
                             promotes_inputs=['initial_states:*', 'final_states:*'])
 
     def configure_ode(self, phase):
@@ -337,7 +336,7 @@ def configure_timeseries_outputs(self, phase):
             The phase object to which this transcription instance applies.
         """
         for timeseries_name, timeseries_options in phase._timeseries.items():
-            timeseries_comp = phase._get_subsystem(f'{timeseries_name}.timeseries_comp')
+            timeseries_comp = phase._get_subsystem(timeseries_name)
             ts_inputs_to_promote = []
             for input_name, src, src_idxs in timeseries_comp._configure_io(timeseries_options):
                 # If the src was added, promote it if it was a state,
@@ -363,12 +362,6 @@ def setup_timeseries_outputs(self, phase):
         gd = self.grid_data
 
         for name, options in phase._timeseries.items():
-            has_expr = False
-            for _, output_options in options['outputs'].items():
-                if output_options['is_expr']:
-                    has_expr = True
-                    break
-
             if options['transcription'] is None:
                 ogd = None
             else:
@@ -378,8 +371,7 @@ def setup_timeseries_outputs(self, phase):
                                                    output_grid_data=ogd,
                                                    output_subset=options['subset'],
                                                    time_units=phase.time_options['units'])
-            timeseries_group = TimeseriesOutputGroup(has_expr=has_expr, timeseries_output_comp=timeseries_comp)
-            phase.add_subsystem(name, subsys=timeseries_group)
+            phase.add_subsystem(name, subsys=timeseries_comp)
 
             phase.connect('dt_dstau', f'{name}.dt_dstau', flat_src_indices=True)
 
@@ -409,7 +401,7 @@ def _get_constraint_kwargs(self, constraint_type, options, phase):
         con_name = constraint_kwargs.pop('constraint_name')
 
         # Determine the path to the variable which we will be constraining
-        var = con_name if options['is_expr'] else options['name']
+        var = options['name']
         var_type = phase.classify_var(var)
 
         # These are the flat indices at a single point in time used
@@ -456,7 +448,6 @@ def _get_constraint_kwargs(self, constraint_type, options, phase):
         con_path = constraint_kwargs.pop('constraint_path')
         constraint_kwargs.pop('shape')
         constraint_kwargs['flat_indices'] = True
-        constraint_kwargs.pop('is_expr')
 
         return con_path, constraint_kwargs
 
diff --git a/dymos/transcriptions/pseudospectral/components/birkhoff_boundary_group.py b/dymos/transcriptions/pseudospectral/components/birkhoff_boundary_group.py
index 9d3138fbb8..2d47f26567 100644
--- a/dymos/transcriptions/pseudospectral/components/birkhoff_boundary_group.py
+++ b/dymos/transcriptions/pseudospectral/components/birkhoff_boundary_group.py
@@ -2,6 +2,7 @@
 import openmdao.api as om
 
 from ...grid_data import GridData
+from ....utils.ode_utils import _make_ode_system
 from dymos._options import options as dymos_options
 
 
@@ -113,12 +114,10 @@ def initialize(self):
                              recordable=False)
         self.options.declare('ode_init_kwargs', types=dict, default={}, recordable=False,
                              desc='Keyword arguments provided when initializing the ODE System')
-        self.options.declare('initial_boundary_constraints', types=list, recordable=False,
-                             desc='Initial boundary constraints from the containing phase.')
-        self.options.declare('final_boundary_constraints', types=list, recordable=False,
-                             desc='Final boundary constraints from the containing phase.')
-        self.options.declare('objectives', types=dict, recordable=False,
-                             desc='Objectives from the containing phase.')
+        self.options.declare('calc_exprs', types=dict, default={},
+                             desc='ODE Expressions from the Phase')
+        self.options.declare('parameter_options', types=dict, default={},
+                             desc='parameter options inherited from the phase')
 
     def setup(self):
         """
@@ -126,19 +125,18 @@ def setup(self):
         """
         ode_class = self.options['ode_class']
         ode_init_kwargs = self.options['ode_init_kwargs']
-        ibcs = self.options['initial_boundary_constraints']
-        fbcs = self.options['final_boundary_constraints']
-        objs = [meta for meta in self.options['objectives'].values()]
 
         self.add_subsystem('boundary_mux', subsys=BirkhoffBoundaryMuxComp(),
                            promotes_inputs=['*'], promotes_outputs=['*'])
 
-        self.add_subsystem('boundary_ode', subsys=ode_class(num_nodes=2, **ode_init_kwargs),
-                           promotes_inputs=['*'], promotes_outputs=['*'])
+        ode = _make_ode_system(ode_class=ode_class,
+                               num_nodes=2,
+                               ode_init_kwargs=ode_init_kwargs,
+                               calc_exprs=self.options['calc_exprs'],
+                               parameter_options=self.options['parameter_options'])
 
-        if any([response['is_expr'] for response in ibcs + fbcs + objs]):
-            self.add_subsystem('boundary_constraint_exec_comp', subsys=om.ExecComp(),
-                               promotes_inputs=['*'], promotes_outputs=['*'])
+        self.add_subsystem('boundary_ode', ode,
+                           promotes_inputs=['*'], promotes_outputs=['*'])
 
     def configure_io(self, phase):
         """
diff --git a/dymos/transcriptions/pseudospectral/components/birkhoff_iter_group.py b/dymos/transcriptions/pseudospectral/components/birkhoff_iter_group.py
index c4a5ab2b4b..149a6e8920 100644
--- a/dymos/transcriptions/pseudospectral/components/birkhoff_iter_group.py
+++ b/dymos/transcriptions/pseudospectral/components/birkhoff_iter_group.py
@@ -5,6 +5,7 @@
 
 from ...grid_data import GridData
 from ....phase.options import TimeOptionsDictionary
+from ....utils.ode_utils import _make_ode_system
 
 
 class BirkhoffIterGroup(om.Group):
@@ -31,12 +32,16 @@ def initialize(self):
                              desc='Dictionary of options for the states.')
         self.options.declare('time_options', types=TimeOptionsDictionary,
                              desc='Options for time in the phase.')
+        self.options.declare('parameter_options', types=dict, default={},
+                             desc='parameter options inherited from the phase')
         self.options.declare('grid_data', types=GridData, desc='Container object for grid info.')
         self.options.declare('ode_class', default=None,
                              desc='Callable that instantiates the ODE system.',
                              recordable=False)
         self.options.declare('ode_init_kwargs', types=dict, default={},
                              desc='Keyword arguments provided when initializing the ODE System')
+        self.options.declare('calc_exprs', types=dict, default={},
+                             desc='ODE Expresions from the Phase')
 
     def setup(self):
         """
@@ -49,7 +54,13 @@ def setup(self):
         ode_class = self.options['ode_class']
         ode_init_kwargs = self.options['ode_init_kwargs']
 
-        self.add_subsystem('ode_all', subsys=ode_class(num_nodes=nn, **ode_init_kwargs))
+        ode = _make_ode_system(ode_class=ode_class,
+                               num_nodes=nn,
+                               ode_init_kwargs=ode_init_kwargs,
+                               calc_exprs=self.options['calc_exprs'],
+                               parameter_options=self.options['parameter_options'])
+
+        self.add_subsystem('ode_all', subsys=ode)
 
         self.add_subsystem('collocation_comp',
                            subsys=BirkhoffDefectComp(grid_data=gd,
diff --git a/dymos/transcriptions/pseudospectral/gauss_lobatto.py b/dymos/transcriptions/pseudospectral/gauss_lobatto.py
index 6389a50f62..3fe0f34fea 100644
--- a/dymos/transcriptions/pseudospectral/gauss_lobatto.py
+++ b/dymos/transcriptions/pseudospectral/gauss_lobatto.py
@@ -8,6 +8,7 @@
 from ...utils.misc import get_rate_units
 from ...utils.introspection import get_promoted_vars, get_targets, get_source_metadata
 from ...utils.indexing import get_src_indices_by_row
+from ...utils.ode_utils import _make_ode_system
 from ..grid_data import GridData, make_subset_map
 
 
@@ -186,9 +187,16 @@ def setup_ode(self, phase):
         grid_data = self.grid_data
         ode_class = phase.options['ode_class']
 
-        kwargs = phase.options['ode_init_kwargs']
-        rhs_disc = ode_class(num_nodes=grid_data.subset_num_nodes['state_disc'], **kwargs)
-        rhs_col = ode_class(num_nodes=grid_data.subset_num_nodes['col'], **kwargs)
+        rhs_disc = _make_ode_system(ode_class=ode_class,
+                                    num_nodes=grid_data.subset_num_nodes['state_disc'],
+                                    ode_init_kwargs=phase.options['ode_init_kwargs'],
+                                    calc_exprs=phase._calc_exprs,
+                                    parameter_options=phase.parameter_options)
+        rhs_col = _make_ode_system(ode_class=ode_class,
+                                   num_nodes=grid_data.subset_num_nodes['col'],
+                                   ode_init_kwargs=phase.options['ode_init_kwargs'],
+                                   calc_exprs=phase._calc_exprs,
+                                   parameter_options=phase.parameter_options)
 
         phase.add_subsystem('rhs_disc', rhs_disc)
 
@@ -297,7 +305,7 @@ def configure_interleave_comp(self, phase):
 
         for timeseries_name, timeseries_options in phase._timeseries.items():
             for ts_output_name, ts_output in timeseries_options['outputs'].items():
-                name = ts_output['name'] if not ts_output['is_expr'] else ts_output_name
+                name = ts_output['name']
                 var_type = phase.classify_var(name)
                 if var_type == 'ode':
                     units = ts_output['units']
diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
index 9d614c760d..2eb94630ab 100644
--- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py
+++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
@@ -2,7 +2,7 @@
 
 import openmdao.api as om
 from ..transcription_base import TranscriptionBase
-from ..common import TimeComp, TimeseriesOutputGroup
+from ..common import TimeComp
 from .components import StateIndependentsComp, StateInterpComp, CollocationComp
 from ..common.timeseries_output_comp import TimeseriesOutputComp
 from ...utils.misc import CoerceDesvar, get_rate_units, reshape_val
@@ -528,12 +528,6 @@ def setup_timeseries_outputs(self, phase):
         gd = self.grid_data
 
         for name, options in phase._timeseries.items():
-            has_expr = False
-            for _, output_options in options['outputs'].items():
-                if output_options['is_expr']:
-                    has_expr = True
-                    break
-
             if options['transcription'] is None:
                 ogd = None
             else:
@@ -543,8 +537,7 @@ def setup_timeseries_outputs(self, phase):
                                                    output_grid_data=ogd,
                                                    output_subset=options['subset'],
                                                    time_units=phase.time_options['units'])
-            timeseries_group = TimeseriesOutputGroup(has_expr=has_expr, timeseries_output_comp=timeseries_comp)
-            phase.add_subsystem(name, subsys=timeseries_group)
+            phase.add_subsystem(name, subsys=timeseries_comp)
 
             phase.connect('dt_dstau', f'{name}.dt_dstau', flat_src_indices=True)
 
diff --git a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
index 4457440e7e..bce1ab4434 100644
--- a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
+++ b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
@@ -6,6 +6,7 @@
 from ...utils.misc import get_rate_units
 from ...utils.introspection import get_promoted_vars, get_targets, get_source_metadata
 from ...utils.indexing import get_src_indices_by_row
+from ...utils.ode_utils import _make_ode_system
 from ..grid_data import GridData
 
 
@@ -119,13 +120,16 @@ def setup_ode(self, phase):
         """
         super(Radau, self).setup_ode(phase)
 
-        ODEClass = phase.options['ode_class']
         grid_data = self.grid_data
 
-        kwargs = phase.options['ode_init_kwargs']
+        ode_sys = _make_ode_system(ode_class=phase.options['ode_class'],
+                                   num_nodes=grid_data.subset_num_nodes['all'],
+                                   ode_init_kwargs=phase.options['ode_init_kwargs'],
+                                   calc_exprs=phase._calc_exprs,
+                                   parameter_options=phase.parameter_options)
+
         phase.add_subsystem('rhs_all',
-                            subsys=ODEClass(num_nodes=grid_data.subset_num_nodes['all'],
-                                            **kwargs))
+                            subsys=ode_sys)
 
     def configure_ode(self, phase):
         """
diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index 5c0c072e8c..6dc1a2017a 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -411,7 +411,7 @@ def configure_timeseries_outputs(self, phase):
             The phase object to which this transcription instance applies.
         """
         for timeseries_name, timeseries_options in phase._timeseries.items():
-            timeseries_comp = phase._get_subsystem(f'{timeseries_name}.timeseries_comp')
+            timeseries_comp = phase._get_subsystem(timeseries_name)
 
             for input_name, src, src_idxs in timeseries_comp._configure_io(timeseries_options):
                 phase.connect(src_name=src,
@@ -465,7 +465,7 @@ def _get_constraint_kwargs(self, constraint_type, options, phase):
         con_name = constraint_kwargs.pop('constraint_name')
 
         # Determine the path to the variable which we will be constraining
-        var = con_name if options['is_expr'] else options['name']
+        var = options['name']
         var_type = phase.classify_var(var)
 
         # These are the flat indices at a single point in time used
@@ -520,7 +520,6 @@ def _get_constraint_kwargs(self, constraint_type, options, phase):
         con_path = constraint_kwargs.pop('constraint_path')
         constraint_kwargs.pop('shape')
         constraint_kwargs['flat_indices'] = True
-        constraint_kwargs.pop('is_expr')
 
         return con_path, constraint_kwargs
 
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index f34ceb9026..45d9de958a 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -1,6 +1,5 @@
 import fnmatch
 from numbers import Number
-import re
 
 import openmdao.api as om
 import numpy as np
@@ -68,12 +67,6 @@ def classify_var(var, time_options, state_options, parameter_options, control_op
     elif var.endswith('_rate2'):
         if var[:-6] in control_options:
             return 'control_rate2'
-    elif timeseries_options is not None:
-        for timeseries in timeseries_options:
-            if var in timeseries_options[timeseries]['outputs']:
-                if timeseries_options[timeseries]['outputs'][var]['is_expr']:
-                    return 'timeseries_exec_comp_output'
-
     return 'ode'
 
 
@@ -171,7 +164,7 @@ def _configure_constraint_introspection(phase):
 
         for con in constraints:
             # Determine the path to the variable which we will be constraining
-            var = con['constraint_name'] if con['is_expr'] else con['name']
+            var = con['name']
             var_type = phase.classify_var(var)
 
             if var != con['constraint_name'] is not None and var_type != 'ode':
@@ -781,16 +774,14 @@ def configure_timeseries_output_introspection(phase):
         not_found = set()
 
         for output_name, output_options in ts_opts['outputs'].items():
-            if output_options['is_expr']:
-                output_meta = phase.timeseries_ec_vars[ts_name][output_name]['meta_data']
-            else:
-                try:
-                    output_meta = transcription._get_timeseries_var_source(output_options['name'],
-                                                                           output_options['output_name'],
-                                                                           phase=phase)
-                except ValueError:
-                    not_found.add(output_name)
-                    continue
+
+            try:
+                output_meta = transcription._get_timeseries_var_source(output_options['name'],
+                                                                       output_options['output_name'],
+                                                                       phase=phase)
+            except ValueError:
+                not_found.add(output_name)
+                raise
 
             output_options['src'] = output_meta['src']
             output_options['src_idxs'] = output_meta['src_idxs']
@@ -810,121 +801,6 @@ def configure_timeseries_output_introspection(phase):
             ts_opts['outputs'].pop(s)
 
 
-def configure_timeseries_expr_introspection(phase):
-    """
-    Introspect the timeseries expressions, add expressions to the exec comp, and make connections to the variables.
-
-    Parameters
-    ----------
-    phase : Phase
-        The phase object whose timeseries outputs are to be introspected.
-    """
-    transcription = phase.options['transcription']
-    num_output_nodes = transcription._get_num_timeseries_nodes()
-    var_names_regex = re.compile(r'([_a-zA-Z]\w*[ ]*\(?:?[.]?)')
-
-    for ts_name, ts_opts in phase._timeseries.items():
-        timeseries_ec = phase._get_subsystem(f'{ts_name}.timeseries_exec_comp')
-        phase.timeseries_ec_vars[ts_name] = {}
-
-        for output_name, output_options in ts_opts['outputs'].items():
-            if output_options['is_expr']:
-
-                expr = output_options['name']
-                expr_lhs = expr.split('=')[0].strip()
-                if '.' in expr_lhs or ':' in expr_lhs:
-                    raise ValueError(
-                        f'{expr_lhs} is an invalid name for the timeseries LHS. Must use a valid variable name')
-                expr_reduced = expr
-
-                units = output_options['units'] if output_options['units'] is not _unspecified else None
-                shape = output_options['shape'] if output_options['shape'] not in _none_or_unspecified else (1,)
-
-                abs_names = [x.strip() for x in re.findall(var_names_regex, expr)
-                             if not x.endswith('(') and not x.endswith(':')]
-                for name in abs_names:
-                    if name.endswith('.'):
-                        idx = abs_names.index(name)
-                        abs_names[idx + 1] = name + abs_names[idx + 1]
-                        abs_names.remove(name)
-                for name in abs_names:
-                    var_name = name.split('.')[-1]
-                    if var_name not in phase.timeseries_ec_vars[ts_name].keys():
-                        phase.timeseries_ec_vars[ts_name][var_name] = {'added_to_ec': False, 'abs_name': name}
-                expr_vars = {}
-                expr_kwargs = output_options['expr_kwargs']
-                kwargs_rel_shape = {}
-                meta = {'units': units, 'shape': shape}
-
-                for var in phase.timeseries_ec_vars[ts_name]:
-                    if var in ('output_name', 'units', 'shape', 'timeseries'):
-                        raise RuntimeError(f'{phase.pathname}: Timeseries expression `{expr}` uses the following'
-                                           f' disallowed variable name in its definition: \'{var}\'.')
-
-                    if phase.timeseries_ec_vars[ts_name][var]['added_to_ec']:
-                        continue
-                    elif var == expr_lhs:
-                        var_units = units
-                        var_shape = shape
-                        meta['src'] = f'{ts_name}.timeseries_exec_comp.{var}'
-                        phase.timeseries_ec_vars[ts_name][var]['added_to_ec'] = True
-                    else:
-                        try:
-                            expr_vars[var] = transcription._get_timeseries_var_source(var, output_name=var, phase=phase)
-                        except ValueError:
-                            expr_vars[var] = transcription._get_timeseries_var_source(
-                                phase.timeseries_ec_vars[ts_name][var]['abs_name'],
-                                output_name=var, phase=phase)
-
-                        var_units = expr_vars[var]['units']
-                        var_shape = expr_vars[var]['shape']
-
-                    if var not in expr_kwargs:
-                        expr_kwargs[var] = {}
-
-                    if 'units' not in expr_kwargs[var]:
-                        expr_kwargs[var]['units'] = var_units
-                    elif var == expr_lhs:
-                        meta['units'] = expr_kwargs[var]['units']
-
-                    if 'shape' in expr_kwargs[var] and var != expr_lhs:
-                        om.issue_warning(f'{phase.pathname}: User-provided shape for timeseries expression input'
-                                         f' `{var}` ignored.\nUsing automatically determined shape'
-                                         f' {(num_output_nodes,) + var_shape}.',
-                                         category=om.UnusedOptionWarning)
-                    elif 'shape' in expr_kwargs[var] and var == expr_lhs:
-                        if output_options['shape'] is not _unspecified:
-                            om.issue_warning(f'{phase.pathname}: User-provided shape for timeseries expression output'
-                                             f' `{var}` using both the\n`shape` keyword argument and the `{var}`'
-                                             f' keyword argument dictionary.\nThe latter will take precedence.',
-                                             category=om.SetupWarning)
-                        var_shape = expr_kwargs[var]['shape']
-                        meta['shape'] = var_shape
-
-                    expr_kwargs[var]['shape'] = var_shape
-                    kwargs_rel_shape[var] = {'units': expr_kwargs[var]['units'],
-                                             'shape': (num_output_nodes,) + expr_kwargs[var]['shape']}
-
-                for input_var in expr_vars:
-                    abs_name = phase.timeseries_ec_vars[ts_name][input_var]['abs_name']
-                    if '.' in abs_name:
-                        expr_reduced = expr_reduced.replace(abs_name, input_var)
-
-                timeseries_ec.add_expr(expr_reduced, **kwargs_rel_shape)
-
-                for var, options in expr_vars.items():
-                    if phase.timeseries_ec_vars[ts_name][var]['added_to_ec']:
-                        continue
-                    else:
-                        phase.connect(src_name=options['src'], tgt_name=f'{ts_name}.timeseries_exec_comp.{var}',
-                                      src_indices=options['src_idxs'])
-
-                    phase.timeseries_ec_vars[ts_name][var]['added_to_ec'] = True
-
-                meta['src_idxs'] = None
-                phase.timeseries_ec_vars[ts_name][expr_lhs]['meta_data'] = meta
-
-
 def filter_outputs(patterns, sys):
     """
     Find all outputs of the given system that match one or more of the strings given in patterns.
@@ -974,7 +850,7 @@ def configure_duration_balance_introspection(phase):
     options = phase.time_options['t_duration_balance_options']
 
     # Determine the path to the variable which we will be constraining
-    var = options['balance_name'] if options['is_expr'] else options['name']
+    var = options['name']
     var_type = phase.classify_var(var)
 
     if var != options['balance_name'] is not None and var_type != 'ode':
diff --git a/dymos/utils/ode_utils.py b/dymos/utils/ode_utils.py
new file mode 100644
index 0000000000..e0dd94baf3
--- /dev/null
+++ b/dymos/utils/ode_utils.py
@@ -0,0 +1,275 @@
+"""
+Utilities for dealing with ODE systems.
+"""
+from copy import deepcopy
+import re
+
+import numpy as np
+
+import openmdao.api as om
+
+from dymos.utils.misc import _unspecified
+
+
+# This regex finds variables and any indices that follow them.
+# 1 - leading underscore or letter (one)
+# 2 - alphanumeric, underscore, period, or colon (zero or more)
+# 3 - end of word is not followed by opening parentheses
+# 4 - optional index specification in brackets
+#                       1111111112222222333333334444444444444444
+var_rgx = re.compile(r'([_a-zA-Z][\w.:]*\b(?!\()(\[[\d:.,-]+\])*)')
+
+
+def _parse_index_string(index_str):
+    """
+    Convert a string representation of array indices to a tuple of slice objects,
+    integer indices, or ellipsis that can be used for array access.
+
+    Args:
+        index_str (str): String representation of indices, e.g. '[1, 3]', '[1:5, 2]', '[..., 0]'
+
+    Returns:
+        tuple: A tuple containing slice objects, integers, or Ellipsis for indexing
+    """
+    if index_str is None:
+        return None
+    # Remove the outer brackets and whitespace
+    clean_str = index_str.strip().strip('[]')
+
+    if not clean_str:
+        return tuple()
+
+    # Split by comma to get each dimension's index
+    dimensions = clean_str.split(',')
+
+    result = []
+    for dim in dimensions:
+        dim = dim.strip()
+
+        # Handle ellipsis
+        if dim == '...' or dim == 'Ellipsis':
+            result.append(Ellipsis)
+            continue
+
+        # Handle full slice ':'
+        if dim == ':':
+            result.append(slice(None))
+            continue
+
+        # Check if it's a slice (contains ':')
+        if ':' in dim:
+            # Handle slice notation (start:stop:step)
+            slice_parts = dim.split(':')
+
+            # Convert each part, handling empty strings
+            start = int(slice_parts[0]) if slice_parts[0].strip() else None
+
+            if len(slice_parts) > 1:
+                stop = int(slice_parts[1]) if slice_parts[1].strip() else None
+            else:
+                stop = None
+
+            if len(slice_parts) > 2:
+                step = int(slice_parts[2]) if slice_parts[2].strip() else None
+            else:
+                step = None
+
+            result.append(slice(start, stop, step))
+        else:
+            # It's a simple index
+            try:
+                result.append(int(dim))
+            except ValueError:
+                # If it's not an integer, it might be a variable
+                result.append(dim)
+
+    return tuple(result)
+
+
+class ExprParser():
+    """
+    A special parser for Expressions in dymos.
+
+    When invoked with a full_path as given by _extract_openmdao_variables,
+    return a legal name for the variable in the expression and an indexer
+    if indices were specified.
+    """
+    def __init__(self):
+        self._unique_int = 0
+
+    def parse(self, full_path, idx_str):
+        """
+        Extract the legal name and any indexer given by the full path representation.
+
+        Parameters
+        ----------
+        full_path : str
+            The full path of the expression variable. This may include
+            dots, colons, and an index specification on the end.
+        idx_str : str
+            The string which follows the path specification in the expression.
+            This begins and ends with brackets [].
+
+        Returns
+        -------
+        tuple(str, tuple or None)
+            The string containing the OpenMDAO-legal name for the variable
+            in the ExecComp expression, as well as any indices specified.
+        """
+        is_ode_rel = False
+        if '.' in full_path:
+            # the path is ODE relative
+            last_path = full_path.split('.')[-1]
+            is_ode_rel = True
+        else:
+            # its a top-level variable or phase variable
+            last_path = full_path
+        legal_name = last_path.split(':')[-1]
+        if is_ode_rel:
+            legal_name += str(self._unique_int)
+            self._unique_int += 1
+        return legal_name, _parse_index_string(idx_str)
+
+
+class ODEGroup(om.Group):
+    """
+    Initialize an ODEGroup to wrap the user's ODE and an ExecComp.
+
+    Parameters
+    ----------
+    ode_class : class
+        The class of ODE used in this Group.
+    num_nodes : int
+        The number of nodes used in the ODE.
+    ode_init_kwargs : dict or None
+        Initialization arguments for the ODE system.
+    calc_exprs : dict
+        The _calc_exprs dictionary of the owning phase instance.
+    parameter_options : dict
+        The parameter_options dictionary of the owning phase instance.
+    """
+    def __init__(self, ode_class, num_nodes, ode_init_kwargs=None, calc_exprs=None, parameter_options=None):
+        super().__init__()
+        self._ode_class = ode_class
+        self._ode_init_kwargs = ode_init_kwargs or {}
+        self._calc_exprs = calc_exprs or {}
+        self._num_nodes = num_nodes
+        self._parameter_options = parameter_options
+
+    def setup(self):
+        """
+        Set up the ODEGroup.
+
+        Adds the user's ODE class and an ExecComp.
+        """
+        ode_class = self._ode_class
+        ode_init_kwargs = self._ode_init_kwargs
+        num_nodes = self._num_nodes
+
+        ode = ode_class(num_nodes=num_nodes, **ode_init_kwargs)
+
+        self.add_subsystem('user_ode', ode, promotes_inputs=['*'], promotes_outputs=['*'])
+        ec = om.ExecComp()
+        self.add_subsystem('exec_comp', ec, promotes_inputs=['*'], promotes_outputs=['*'])
+
+    def configure(self):
+        """
+        Setup up connections and promotions in the ODEGroup.
+        """
+        num_nodes = self._num_nodes
+
+        seen_kwargs = set()
+        parser = ExprParser()
+        ec = self._get_subsystem('exec_comp')
+
+        for expr, expr_kwargs in self._calc_exprs.items():
+            common_units = _unspecified
+            common_shape = _unspecified
+            if 'units' in expr_kwargs:
+                # units are set throughout this expression
+                common_units = expr_kwargs['units']
+            if 'shape' in expr_kwargs:
+                common_shape = expr_kwargs['shape']
+
+            _expr_kwargs = {}
+            output_var, rhs = [s.strip() for s in expr.split('=')]
+            _expr_kwargs[output_var] = deepcopy(expr_kwargs.get(output_var, {}))
+
+            if 'shape' not in _expr_kwargs[output_var]:
+                if common_shape is not _unspecified:
+                    _expr_kwargs[output_var]['shape'] = (num_nodes,) + common_shape
+                else:
+                    _expr_kwargs[output_var]['shape'] = (num_nodes,)
+            else:
+                # Assume given shape is at a per-node basis
+                _expr_kwargs[output_var]['shape'] = (num_nodes,) + _expr_kwargs[output_var]['shape']
+
+            if 'units' not in _expr_kwargs[output_var]:
+                if common_units is not _unspecified:
+                    _expr_kwargs[output_var]['units'] = common_units
+
+            scalar_src_idxs = []
+            param_options = self._parameter_options or {}
+            scalar_sources = list(param_options.keys()) + ['t_initial', 't_duration', 't_final']
+
+            for rel_path, idx_str in re.findall(var_rgx, rhs):
+                exec_var_name, src_idxs = parser.parse(rel_path, idx_str)
+                expr = expr.replace(rel_path, exec_var_name)
+                if '.' in rel_path:
+                    self.connect(rel_path, exec_var_name, src_indices=src_idxs)
+
+                # Only provide kwargs for things that we havent already done so.
+                if exec_var_name not in seen_kwargs:
+                    # Use deepcopy so we don't accidentally permanently set the shape here when we assign it.
+                    _expr_kwargs[exec_var_name] = deepcopy(expr_kwargs.get(rel_path, {}))
+                    if 'shape' not in _expr_kwargs[exec_var_name]:
+                        if common_shape is not _unspecified:
+                            _expr_kwargs[exec_var_name]['shape'] = (num_nodes,) + common_shape
+                        else:
+                            _expr_kwargs[exec_var_name]['shape'] = (num_nodes,)
+                    else:
+                        _expr_kwargs[exec_var_name]['shape'] = (num_nodes,) + _expr_kwargs[exec_var_name]['shape']
+
+                    if exec_var_name in scalar_sources:
+                        scalar_src_idxs.append(exec_var_name)
+
+            seen_kwargs |= _expr_kwargs.keys()
+            ec.add_expr(expr, **_expr_kwargs)
+            if scalar_src_idxs:
+                self.promotes('exec_comp', inputs=scalar_src_idxs,
+                              src_indices=np.zeros(num_nodes, dtype=int))
+
+
+def _make_ode_system(ode_class, num_nodes, ode_init_kwargs=None, calc_exprs=None, parameter_options=None):
+    """
+    Instantiate the ODE system, optionally including an ExecComp.
+
+    Parameters
+    ----------
+    ode_class : System class
+        The ODE class provided to the dymos phase.
+    num_nodes : int
+        The number of nodes at which the ODE is instantiated.
+    ode_init_kwargs : dict or None
+        Additional arguments we need to instantiate the ODE.
+    calc_exprs : dict
+        A dictionary keyed by exec comp expresions whose associated values
+        are the metadata associated with the expression.
+
+    Returns
+    -------
+    ode : System
+        The instantiation of ode_class.  If `calc_exprs` were given, it
+        returns a group that wraps the instantiated ode_class, as well
+        as an exec comp used to evaluate the expressions.
+    """
+    _kwargs = ode_init_kwargs or {}
+    if not calc_exprs:
+        return ode_class(num_nodes=num_nodes, **_kwargs)
+    else:
+        ode_group = ODEGroup(ode_class,
+                             num_nodes=num_nodes,
+                             ode_init_kwargs=ode_init_kwargs,
+                             calc_exprs=calc_exprs,
+                             parameter_options=parameter_options)
+        return ode_group
diff --git a/dymos/visualization/timeseries_plots.py b/dymos/visualization/timeseries_plots.py
index 37a02b30dc..52af404bd5 100644
--- a/dymos/visualization/timeseries_plots.py
+++ b/dymos/visualization/timeseries_plots.py
@@ -397,19 +397,17 @@ def timeseries_plots(solution_recorder_filename, simulation_record_file=None, pl
                 break
         # get the time units first so they can be associated with all the variables
         for timeseries_node_child in timeseries_node['children']:
-            for timeseries_node_g_child in timeseries_node_child['children']:
-                if timeseries_node_g_child['name'] == 'time':
-                    units_for_time = timeseries_node_g_child['units']
-                    break
+            if timeseries_node_child['name'] == 'time':
+                units_for_time = timeseries_node_child['units']
+                break
         for timeseries_node_child in timeseries_node['children']:
             # plot everything in the timeseries except input_values. Also not time since that
             #   is the independent variable in these plot
-            for timeseries_node_g_child in timeseries_node_child['children']:
-                if not timeseries_node_g_child['name'].startswith("input_values:")\
-                        and not timeseries_node_g_child['name'] == 'time':
-                    varname = timeseries_node_g_child['name']
-                    var_units[varname] = timeseries_node_g_child['units']
-                    time_units[varname] = units_for_time
+            if not timeseries_node_child['name'].startswith("input_values:")\
+                    and not timeseries_node_child['name'] == 'time':
+                varname = timeseries_node_child['name']
+                var_units[varname] = timeseries_node_child['units']
+                time_units[varname] = units_for_time
 
     # Check to see if there is anything to plot
     if len(var_units) == 0:


Test patch:
diff --git a/docs/dymos_book/test/test_jupyter_output_linting.py b/docs/dymos_book/test/test_jupyter_output_linting.py
index 0aa55ea8d2..3956d780e3 100644
--- a/docs/dymos_book/test/test_jupyter_output_linting.py
+++ b/docs/dymos_book/test/test_jupyter_output_linting.py
@@ -38,7 +38,6 @@ class LintJupyterOutputsTestCase(unittest.TestCase):
     def test_output(self):
 
         for file in _get_files():
-            print(file)
             with self.subTest(file):
                 with open(file) as f:
                     json_data = json.load(f)
diff --git a/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py b/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py
index e4ccd25c7f..e9ef1d4002 100644
--- a/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py
+++ b/dymos/examples/balanced_field/doc/test_doc_balanced_field_length.py
@@ -200,9 +200,12 @@ def test_balanced_field_length_for_docs(self):
         rto.add_objective('r', loc='final', ref=1.0)
 
         for phase_name, phase in traj._phases.items():
-            phase.add_timeseries_output('T_nominal', output_name='T')
-            phase.add_timeseries_output('T_engine_out', output_name='T')
-            phase.add_timeseries_output('T_shutdown', output_name='T')
+            if 'T_nominal' in phase.parameter_options:
+                phase.add_timeseries_output('T_nominal', output_name='T')
+            if 'T_engine_out' in phase.parameter_options:
+                phase.add_timeseries_output('T_engine_out', output_name='T')
+            if 'T_shutdown' in phase.parameter_options:
+                phase.add_timeseries_output('T_shutdown', output_name='T')
             phase.add_timeseries_output('alpha')
 
         #
diff --git a/dymos/examples/cannonball/test/test_cannonball_connected_initial_states.py b/dymos/examples/cannonball/test/test_cannonball_connected_initial_states.py
index 9e46256ecc..c64984f651 100644
--- a/dymos/examples/cannonball/test/test_cannonball_connected_initial_states.py
+++ b/dymos/examples/cannonball/test/test_cannonball_connected_initial_states.py
@@ -53,7 +53,7 @@ def test_cannonball_connected_initial_states(self):
         ascent.add_parameter('m', units='kg', static_targets=True)
 
         # Limit the muzzle energy
-        ascent.add_boundary_constraint('ke = 0.5 * m * v**2', loc='initial',
+        ascent.add_boundary_constraint('ke2 = 0.5 * m * v**2', loc='initial',
                                        upper=400000, lower=0, ref=100000)
 
         # Second Phase (descent)
diff --git a/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py b/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py
index fa8f5b497a..6ff69399b3 100644
--- a/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py
+++ b/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py
@@ -52,7 +52,7 @@ def test_cannonball_design_boundary_constraint_expression(self):
         ascent.add_parameter('m', units='kg', dynamic=False)
 
         # Limit the muzzle energy
-        ascent.add_boundary_constraint('ke = 0.5 * m * v**2', loc='initial',
+        ascent.add_boundary_constraint('ke2 = 0.5 * m * v**2', loc='initial',
                                        upper=400000, lower=0, ref=100000)
 
         # Second Phase (descent)
diff --git a/dymos/phase/test/test_add_boundary_constraint.py b/dymos/phase/test/test_add_boundary_constraint.py
index 8df2376084..60b37651a6 100644
--- a/dymos/phase/test/test_add_boundary_constraint.py
+++ b/dymos/phase/test/test_add_boundary_constraint.py
@@ -99,7 +99,6 @@ def test_invalid_expression(self):
 
         phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='s')
 
-        # can't fix final position if you're solving the segments
         phase.add_state('pos',
                         rate_source='pos_dot', units='m',
                         fix_initial=True)
@@ -115,11 +114,11 @@ def test_invalid_expression(self):
         phase.add_parameter('g', units='m/s**2', val=9.80665, opt=False)
 
         # test add_boundary_constraint with arrays:
-        expected = 'The expression provided `pos**2` has invalid format. ' \
-                   'Expression may be a single variable or an equation ' \
-                   'of the form `constraint_name = func(vars)`'
+        phase.add_boundary_constraint(name='pos**2', loc='final', equals=np.array([10, 5]))
+
+        expected = "Unable to find the source 'pos**2' in the ODE."
         with self.assertRaises(ValueError) as e:
-            phase.add_boundary_constraint(name='pos**2', loc='final', equals=np.array([10, 5]))
+            p.setup()
 
         self.assertEqual(expected, str(e.exception))
 
@@ -165,8 +164,8 @@ def test_duplicate_name(self):
         with self.assertRaises(ValueError) as e:
             phase.add_boundary_constraint(name='pos=v**2', loc='final', equals=np.array([10, 5]))
 
-        expected = 'Cannot add new final boundary constraint named `pos` and indices None.' \
-                   ' The name `pos` is already in use as a final boundary constraint'
+        expected = 'Cannot add new final boundary constraint for variable `pos` and indices None.' \
+                   ' One already exists.'
         self.assertEqual(expected, str(e.exception))
 
     def test_duplicate_constraint(self):
diff --git a/dymos/phase/test/test_add_path_constraint.py b/dymos/phase/test/test_add_path_constraint.py
index 345c664ccf..55e763aed0 100644
--- a/dymos/phase/test/test_add_path_constraint.py
+++ b/dymos/phase/test/test_add_path_constraint.py
@@ -85,12 +85,12 @@ def test_invalid_expression(self):
 
         phase.add_parameter('g', units='m/s**2', val=9.80665, opt=False)
 
+        phase.add_path_constraint(name='pos**2', lower=np.array([10, 5]))
+
         # test add_boundary_constraint with arrays:
-        expected = 'The expression provided `pos**2` has invalid format. ' \
-                   'Expression may be a single variable or an equation ' \
-                   'of the form `constraint_name = func(vars)`'
+        expected = "Unable to find the source 'pos**2' in the ODE."
         with self.assertRaises(ValueError) as e:
-            phase.add_path_constraint(name='pos**2', equals=np.array([10, 5]))
+            p.setup()
 
         self.assertEqual(expected, str(e.exception))
 
diff --git a/dymos/phase/test/test_analytic_phase.py b/dymos/phase/test/test_analytic_phase.py
index 06416ddeb4..a2e0580b8f 100644
--- a/dymos/phase/test/test_analytic_phase.py
+++ b/dymos/phase/test/test_analytic_phase.py
@@ -96,7 +96,27 @@ def test_simple_ivp_system(self):
 
         assert_near_equal(y, expected_y)
 
-        p.check_partials(compact_print=True, method='cs')
+        p.check_partials(compact_print=True, method='cs', out_stream=None)
+
+    def test_simple_ivp_system(self):
+
+        p = om.Problem()
+        p.model.add_subsystem('simple_ivp_solution', SimpleIVPSolution(num_nodes=100), promotes=['*'])
+        p.setup(force_alloc_complex=True)
+
+        p.set_val('t', np.linspace(0, 2, 100))
+        p.set_val('y0', 0.5)
+
+        p.run_model()
+
+        t = p.get_val('t')
+        y = p.get_val('y')
+
+        expected_y = t ** 2 + 2 * t + 1 - 0.5 * np.exp(t)
+
+        assert_near_equal(y, expected_y)
+
+        p.check_partials(compact_print=True, method='cs', out_stream=None)
 
     def test_simple_ivp(self):
 
@@ -126,6 +146,38 @@ def test_simple_ivp(self):
 
         assert_near_equal(y, expected)
 
+    def test_simple_ivp_calc_expr(self):
+
+        p = om.Problem()
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+        phase = dm.AnalyticPhase(ode_class=SimpleIVPSolution, num_nodes=10)
+
+        phase.add_calc_expr('y2 = y**2')
+
+        traj.add_phase('phase', phase)
+
+        phase.set_time_options(units='s', targets=['t'], fix_initial=True, fix_duration=True)
+        phase.add_state('y')
+        phase.add_parameter('y0', opt=False)
+
+        p.setup()
+
+        p.set_val('traj.phase.t_initial', 0.0, units='s')
+        p.set_val('traj.phase.t_duration', 2.0, units='s')
+        p.set_val('traj.phase.parameters:y0', 0.5, units='unitless')
+
+        p.run_model()
+
+        t = p.get_val('traj.phase.timeseries.time', units='s')
+        y = p.get_val('traj.phase.timeseries.y', units='unitless')
+        y2 = p.get_val('traj.phase.timeseries.y2')
+
+        expected = t ** 2 + 2 * t + 1 - 0.5 * np.exp(t)
+
+        assert_near_equal(y, expected)
+        assert_near_equal(y2, y*y)
+
     def test_simple_bvp(self):
 
         p = om.Problem()
@@ -299,8 +351,8 @@ def test_timeseries_expr(self):
 
         phase.set_time_options(units='s', targets=['t'], fix_initial=True, fix_duration=True)
         phase.add_state('y')
-        phase.add_parameter('y0', opt=False)
-        phase.add_timeseries_output('z = y0 + y**2')
+        phase.add_parameter('y0', units='unitless', opt=False)
+        phase.add_timeseries_output('z = y0 + y**2', y0={'units': 'unitless'})
 
         p.setup()
 
diff --git a/dymos/phase/test/test_implicit_duration.py b/dymos/phase/test/test_implicit_duration.py
index ba244c1008..046b135250 100644
--- a/dymos/phase/test/test_implicit_duration.py
+++ b/dymos/phase/test/test_implicit_duration.py
@@ -74,6 +74,9 @@ def _make_simple_problem(tx, input_duration=False, fix_duration=False, expr_end=
         p.driver = om.ScipyOptimizeDriver(optimizer='SLSQP')
         p.driver.declare_coloring(tol=1.0E-12)
 
+        phase.nonlinear_solver = om.NewtonSolver(solve_subsystems=True)
+        phase.linear_solver = om.DirectSolver()
+
         p.setup()
 
         p.set_val('traj.phase.t_initial', 0)
@@ -100,6 +103,9 @@ def _make_matrix_problem(tx, shape_error=False):
 
         phase.set_simulate_options(rtol=1.0E-9, atol=1.0E-9)
 
+        phase.nonlinear_solver = om.NewtonSolver(solve_subsystems=True)
+        phase.linear_solver = om.DirectSolver()
+
         p.driver = om.ScipyOptimizeDriver(optimizer='SLSQP')
         p.driver.declare_coloring(tol=1.0E-12)
 
@@ -181,6 +187,7 @@ def test_implicit_duration_matrix_state(self):
 
         assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
 
+    @unittest.skip('ODE Expressions currently cannot be provided in balances.')
     @require_pyoptsparse(optimizer='IPOPT')
     def test_implicit_duration_radau_expr_condition(self):
         tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
@@ -194,6 +201,7 @@ def test_implicit_duration_radau_expr_condition(self):
 
         print((p.get_val('traj.phase.timeseries.v')[-1]) ** 2 / 2)
 
+    @unittest.skip('ODE Expressions currently cannot be provided in balances.')
     def test_implicit_duration_gl_expr_condition(self):
         tx = dm.GaussLobatto(num_segments=12, order=3, solve_segments=False)
 
diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index 4ef4c8bef6..b53bfad8ea 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -385,6 +385,7 @@ def min_time_climb(num_seg=3, transcription_class=dm.Radau, transcription_order=
     return p
 
 
+@use_tempdirs
 class TestDuplicateTimeseriesGlobName(unittest.TestCase):
 
     def test_duplicate_timeseries_glob_name_radau(self):
@@ -406,6 +407,7 @@ def test_duplicate_timeseries_glob_name_gl(self):
         self.assertEqual(str(e.exception), msg)
 
 
+@use_tempdirs
 class TestTimeseriesExprBrachistochrone(unittest.TestCase):
 
     @staticmethod
@@ -444,22 +446,14 @@ def make_problem_brachistochrone(transcription, polynomial_control=False):
 
         p.model.options['assembled_jac_type'] = 'csc'
         p.model.linear_solver = om.DirectSolver()
-        p.setup(check=True)
-
-        p['phase0.t_initial'] = 0.0
-        p['phase0.t_duration'] = 2.0
-
-        if isinstance(transcription, dm.ExplicitShooting):
-            p['phase0.initial_states:x'] = 0
-            p['phase0.initial_states:y'] = 10
-            p['phase0.initial_states:v'] = 0
-        else:
-            p['phase0.states:x'] = phase.interp('x', [0, 10])
-            p['phase0.states:y'] = phase.interp('y', [10, 5])
-            p['phase0.states:v'] = phase.interp('v', [0, 9.9])
-
-        p[f'phase0.{control_name}'] = phase.interp('theta', [5, 100])
-        p['phase0.parameters:g'] = 9.80665
+        p.setup()
+
+        phase.set_time_val(initial=0.0, duration=2.0)
+        phase.set_state_val('x', [0, 10])
+        phase.set_state_val('y', [10, 5])
+        phase.set_state_val('v', [0, 9.9])
+        phase.set_control_val('theta', [5, 100])
+        phase.set_parameter_val('g', 9.80665)
         return p
 
     def test_invalid_expr_var(self):
@@ -484,11 +478,11 @@ def test_invalid_expr_var(self):
 
         phase.add_timeseries_output('k=units', units='m**2')
 
-        with self.assertRaises(RuntimeError) as e:
+        with self.assertRaises(NameError) as e:
             p.setup(check=True)
 
-        expected = "phase0: Timeseries expression `k=units` uses the following disallowed variable " \
-                   "name in its definition: 'units'."
+        expected = ("'phase0.rhs_all.exec_comp' <class ExecComp>: cannot use variable "
+                    "name 'units' because it's a reserved keyword.")
         self.assertEqual(str(e.exception), expected)
 
     def test_input_units(self):
@@ -551,69 +545,26 @@ def test_output_units(self):
 
         phase.add_parameter('g', opt=False, units='m/s**2', val=9.80665, include_timeseries=True)
 
-        phase.add_timeseries_output('sin_theta=sin(theta)', units='unitless', theta={'units': 'rad'})
+        phase.add_timeseries_output('sin_theta=sin(theta)', shape=(1,), units='unitless', theta={'units': 'rad', 'shape': (1,)})
         phase.add_timeseries_output('sin_theta2=sin(theta)', sin_theta2={'units': 'unitless'})
         phase.add_timeseries_output('sin_theta3=sin(theta)', shape=(1,), sin_theta3={'shape': (1,)})
 
-        expected = 'phase0: User-provided shape for timeseries expression output `sin_theta3` using ' \
-                   'both the\n`shape` keyword argument and the `sin_theta3` keyword argument dictionary.\n' \
-                   'The latter will take precedence.'
-
-        with warnings.catch_warnings(record=True) as ctx:
-            warnings.simplefilter('always')
-            p.setup(check=True)
-
-        self.assertIn(expected, [str(warning.message) for warning in ctx])
+        p.setup()
 
-        p['phase0.t_initial'] = 0.0
-        p['phase0.t_duration'] = 2.0
-
-        p['phase0.states:x'] = phase.interp('x', [0, 10])
-        p['phase0.states:y'] = phase.interp('y', [10, 5])
-        p['phase0.states:v'] = phase.interp('v', [0, 9.9])
-        p['phase0.controls:theta'] = phase.interp('theta', [5, 100])
-        p['phase0.parameters:g'] = 9.80665
+        phase.set_time_val(initial=0.0, duration=2.0)
+        phase.set_state_val('x', [0, 10])
+        phase.set_state_val('y', [10, 5])
+        phase.set_state_val('v', [0, 9.9])
+        phase.set_control_val('theta', [5, 100])
+        phase.set_parameter_val('g', 9.80665)
 
         p.run_model()
 
-        from dymos.utils.introspection import get_promoted_vars
-
-        sin_theta_meta = get_promoted_vars(phase, iotypes='output', metadata_keys=('units',))['timeseries.sin_theta']
-        sin_theta2_meta = get_promoted_vars(phase, iotypes='output', metadata_keys=('units',))['timeseries.sin_theta2']
-
-        self.assertEqual('unitless', sin_theta_meta['units'])
-        self.assertEqual('unitless', sin_theta2_meta['units'])
-
-    def test_ignore_shape(self):
-        p = om.Problem(model=om.Group())
-
-        tx = dm.Radau(num_segments=5, order=3, compressed=True)
+        var_data = p.model.list_vars(units=True, shape=True, out_stream=None)
+        var_data = {meta['prom_name']: meta for abs_path, meta in var_data}
 
-        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
-
-        p.model.add_subsystem('phase0', phase)
-
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
-
-        phase.set_state_options('x', fix_initial=True)
-        phase.set_state_options('y', fix_initial=True)
-        phase.set_state_options('v', fix_initial=True)
-
-        phase.add_control('theta', continuity=True, rate_continuity=True, opt=True,
-                          units='deg', lower=0.01, upper=179.9, ref=1, ref0=0)
-
-        phase.add_parameter('g', opt=False, units='m/s**2', val=9.80665, include_timeseries=True)
-
-        phase.add_timeseries_output('sin_theta=sin(theta)', units='unitless', theta={'units': 'rad', 'shape': (1,)})
-
-        expected = 'phase0: User-provided shape for timeseries expression input `theta` ignored.\n' \
-                   'Using automatically determined shape (20, 1).'
-
-        with warnings.catch_warnings(record=True) as ctx:
-            warnings.simplefilter('always')
-            p.setup(check=True)
-
-        self.assertIn(expected, [str(warning.message) for warning in ctx])
+        self.assertEqual('unitless', var_data['phase0.timeseries.sin_theta']['units'])
+        self.assertEqual('unitless', var_data['phase0.timeseries.sin_theta']['units'])
 
     def test_timeseries_expr_radau(self):
         tx = dm.Radau(num_segments=5, order=3, compressed=True)
@@ -666,41 +617,63 @@ def test_timeseries_expr_gl(self):
         assert_near_equal(z_computed, z_ts, tolerance=1e-12)
         assert_near_equal(f_computed, f_ts, tolerance=1e-12)
 
-    def test_timeseries_expr_gl_polynomial_control(self):
-        tx = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
-        p = self.make_problem_brachistochrone(transcription=tx, polynomial_control=True)
+    def test_timeseries_expr_birkhoff(self):
+        tx = dm.Birkhoff(num_nodes=21)
+        p = self.make_problem_brachistochrone(transcription=tx)
         p.run_driver()
-
         x = p.get_val('phase0.timeseries.x')
         y = p.get_val('phase0.timeseries.y')
         theta = p.get_val('phase0.timeseries.theta')
         g = p.get_val('phase0.timeseries.g')
 
         z_computed = x * y + x**2
-        f_computed = 3 * g * np.cos(theta)**2
+        f_computed = 3 * g * np.cos(theta) ** 2
+
         z_ts = p.get_val('phase0.timeseries.z')
         f_ts = p.get_val('phase0.timeseries.f')
         assert_near_equal(z_computed, z_ts, tolerance=1e-12)
         assert_near_equal(f_computed, f_ts, tolerance=1e-12)
 
-    def test_timeseries_expr_explicit_shooting(self):
-        tx = dm.ExplicitShooting(grid=dm.GaussLobattoGrid(num_segments=3, nodes_per_seg=5, compressed=True))
-
-        p = self.make_problem_brachistochrone(transcription=tx)
+    def test_timeseries_expr_gl_polynomial_control(self):
+        tx = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
+        p = self.make_problem_brachistochrone(transcription=tx, polynomial_control=True)
         p.run_driver()
+
         x = p.get_val('phase0.timeseries.x')
         y = p.get_val('phase0.timeseries.y')
         theta = p.get_val('phase0.timeseries.theta')
         g = p.get_val('phase0.timeseries.g')
 
         z_computed = x * y + x**2
-        f_computed = 3 * g * np.cos(theta) ** 2
-
+        f_computed = 3 * g * np.cos(theta)**2
         z_ts = p.get_val('phase0.timeseries.z')
         f_ts = p.get_val('phase0.timeseries.f')
         assert_near_equal(z_computed, z_ts, tolerance=1e-12)
         assert_near_equal(f_computed, f_ts, tolerance=1e-12)
 
+    def test_timeseries_expr_explicit_shooting(self):
+
+        import warnings
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", category=om.UnitsWarning)
+            tx = dm.ExplicitShooting(grid=dm.GaussLobattoGrid(num_segments=3, nodes_per_seg=5, compressed=True))
+
+            p = self.make_problem_brachistochrone(transcription=tx)
+            p.run_driver()
+            x = p.get_val('phase0.timeseries.x')
+            y = p.get_val('phase0.timeseries.y')
+            theta = p.get_val('phase0.timeseries.theta')
+            g = p.get_val('phase0.timeseries.g')
+
+            z_computed = x * y + x**2
+            f_computed = 3 * g * np.cos(theta) ** 2
+
+            z_ts = p.get_val('phase0.timeseries.z')
+            f_ts = p.get_val('phase0.timeseries.f')
+            assert_near_equal(z_computed, z_ts, tolerance=1e-12)
+            assert_near_equal(f_computed, f_ts, tolerance=1e-12)
+
     def test_timeseries_expr_explicit_shooting_polynomial_controls(self):
         tx = dm.ExplicitShooting(grid=dm.GaussLobattoGrid(num_segments=3, nodes_per_seg=5, compressed=True))
 
@@ -754,7 +727,8 @@ def test_timeseries_expr_solve_ivp_polynomial_controls(self):
         assert_near_equal(f_ts[-1], f_sim[-1], tolerance=1e-3)
 
 
-class TestTimeseriesMinTimeClimb(unittest.TestCase):
+@use_tempdirs
+class TestTimeseriesExprMinTimeClimb(unittest.TestCase):
 
     def test_timeseries_expr_radau(self):
         p = min_time_climb(transcription_class=dm.Radau, num_seg=12, transcription_order=3, timeseries_expr=True)
diff --git a/dymos/test/test_load_case.py b/dymos/test/test_load_case.py
index b336ca87ac..4f663c7e0f 100644
--- a/dymos/test/test_load_case.py
+++ b/dymos/test/test_load_case.py
@@ -153,13 +153,13 @@ def test_load_case_lgl_to_radau(self):
         outputs = dict([(o[0], o[1]) for o in case.list_outputs(units=True, shape=True,
                                                                 out_stream=None)])
 
-        time_val = outputs['phase0.timeseries.timeseries_comp.time']['val']
-        theta_val = outputs['phase0.timeseries.timeseries_comp.theta']['val']
+        time_val = outputs['phase0.timeseries.time']['val']
+        theta_val = outputs['phase0.timeseries.theta']['val']
 
         time_val_uniq, idx = np.unique(time_val, return_index=True)
         theta_val_uniq = theta_val[idx]
 
-        assert_near_equal(q['phase0.timeseries.timeseries_comp.theta'],
+        assert_near_equal(q['phase0.timeseries.theta'],
                           q.model.phase0.interp(xs=time_val_uniq, ys=theta_val_uniq, nodes='all'),
                           tolerance=1.0E-3)
 
diff --git a/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py b/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py
index fad5d62c3a..0a729c3db5 100644
--- a/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py
@@ -544,14 +544,12 @@ def test_explicit_shooting_unknown_timeseries(self):
 
             phase.add_objective('time', loc='final')
 
-            msg = "phase0: The following timeseries outputs were requested but not found in the " \
-                  "ODE: foo"
+            msg = "Unable to find the source 'foo' in the ODE."
 
-            with warnings.catch_warnings(record=True) as ctx:
-                warnings.simplefilter('always')
-                prob.setup(force_alloc_complex=True)
+            with self.assertRaises(ValueError) as e:
+                prob.setup()
 
-            self.assertIn(msg, [str(w.message) for w in ctx])
+            self.assertIn(msg, str(e.exception))
 
     def test_brachistochrone_static_gravity_explicit_shooting(self):
         import openmdao.api as om
diff --git a/dymos/transcriptions/pseudospectral/test/test_ps_base.py b/dymos/transcriptions/pseudospectral/test/test_ps_base.py
index 145adf61b5..47791a5721 100644
--- a/dymos/transcriptions/pseudospectral/test/test_ps_base.py
+++ b/dymos/transcriptions/pseudospectral/test/test_ps_base.py
@@ -87,7 +87,7 @@ def test_add_state_units_none(self):
         p = make_problem(transcription=GaussLobatto, num_segments=10, order=3, compressed=True)
         p.run_model()
 
-        io_meta = p.model.traj.phases.phase.timeseries.timeseries_comp.get_io_metadata(
+        io_meta = p.model.traj.phases.phase.timeseries.get_io_metadata(
             iotypes=('output',), get_remote=True)
         self.assertEqual(io_meta['x']['units'], None)