diff --git a/.github/workflows/dymos_tests_workflow.yml b/.github/workflows/dymos_tests_workflow.yml
index 9c4fcff46..edf3e4c98 100644
--- a/.github/workflows/dymos_tests_workflow.yml
+++ b/.github/workflows/dymos_tests_workflow.yml
@@ -285,7 +285,7 @@ jobs:
             if [[ "${{ matrix.SNOPT }}" ]]; then
               echo "SNOPT ${{ matrix.SNOPT }} was requested but is not available on conda-forge"
             fi
-
+            conda install -c conda-forge lapack
             conda install -c conda-forge pyoptsparse
           else
             pip install git+https://github.com/OpenMDAO/build_pyoptsparse
@@ -441,9 +441,9 @@ jobs:
           testflo -b benchmark --pre_announce
           cd $HOME
           if [[ "${{ matrix.NAME }}" != "latest" ]]; then
-            testflo dymos -n 4 --pre_announce --show_skipped --durations 20 --coverage --coverpkg dymos
+            testflo dymos -n 1 --pre_announce --show_skipped --durations 20 --coverage --coverpkg dymos
           else
-            testflo dymos -n 4 --pre_announce --show_skipped --durations 20
+            testflo dymos -n 1 --pre_announce --show_skipped --durations 20
           fi
       - name: Submit coverage
         if: github.event_name != 'workflow_dispatch'
diff --git a/benchmark/benchmark_racecar.py b/benchmark/benchmark_racecar.py
index 413c7c720..411725c97 100644
--- a/benchmark/benchmark_racecar.py
+++ b/benchmark/benchmark_racecar.py
@@ -140,22 +140,22 @@ def _run_racecar_problem(transcription, timeseries=False, make_plots=False):
     # Now that the OpenMDAO problem is setup, we can set the values of the states.
 
     # States
-    # non-zero velocity in order to protect against 1/0 errors.
-    p.set_val('traj.phase0.states:V', phase.interp('V', [20, 20]), units='m/s')
-    p.set_val('traj.phase0.states:lambda', phase.interp('lambda', [0.0, 0.0]), units='rad')
-    # all other states start at 0
-    p.set_val('traj.phase0.states:omega', phase.interp('omega', [0.0, 0.0]), units='rad/s')
-    p.set_val('traj.phase0.states:alpha', phase.interp('alpha', [0.0, 0.0]), units='rad')
-    p.set_val('traj.phase0.states:ax', phase.interp('ax', [0.0, 0.0]), units='m/s**2')
-    p.set_val('traj.phase0.states:ay', phase.interp('ay', [0.0, 0.0]), units='m/s**2')
-    p.set_val('traj.phase0.states:n', phase.interp('n', [0.0, 0.0]), units='m')
+    # Nonzero velocity to avoid division by zero errors
+    phase.set_state_val('V', 20.0, units='m/s')
+    # All other states start at 0
+    phase.set_state_val('lambda', 0.0, units='rad')
+    phase.set_state_val('omega', 0.0, units='rad/s')
+    phase.set_state_val('alpha', 0.0, units='rad')
+    phase.set_state_val('ax', 0.0, units='m/s**2')
+    phase.set_state_val('ay', 0.0, units='m/s**2')
+    phase.set_state_val('n', 0.0, units='m')
     # initial guess for what the final time should be
-    p.set_val('traj.phase0.states:t', phase.interp('t', [0.0, 100.0]), units='s')
+    phase.set_state_val('t', [0.0, 100.0], units='s')
 
     # Controls
-    p.set_val('traj.phase0.controls:delta', phase.interp('delta', [0.0, 0.0]), units='rad')
-    p.set_val('traj.phase0.controls:thrust', phase.interp('thrust', [0.1, 0.1]), units=None)
     # a small amount of thrust can speed up convergence
+    phase.set_control_val('delta', 0.0, units='rad')
+    phase.set_control_val('thrust', 0.1, units=None)
 
     dm.run_problem(p, run_driver=True, simulate=False, make_plots=make_plots)
     print('Optimization finished')
@@ -183,4 +183,4 @@ def benchmark_radau_timeseries(self):
 
 
 if __name__ == '__main__':
-    _run_racecar_problem(dm.GaussLobatto, timeseries=False, make_plots=True)
+    _run_racecar_problem(dm.Radau, timeseries=False, make_plots=True)
diff --git a/dymos/__init__.py b/dymos/__init__.py
index cbb6bb70c..36c499204 100644
--- a/dymos/__init__.py
+++ b/dymos/__init__.py
@@ -2,7 +2,7 @@
 
 
 from .phase import Phase, AnalyticPhase
-from .transcriptions import GaussLobatto, Radau, ExplicitShooting, Analytic, Birkhoff
+from .transcriptions import GaussLobatto, Radau, ExplicitShooting, Analytic, Birkhoff, RadauDeprecated
 from .transcriptions.grid_data import GaussLobattoGrid, RadauGrid, UniformGrid, BirkhoffGrid
 from .trajectory.trajectory import Trajectory
 from .run_problem import run_problem
diff --git a/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py b/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py
index a5337817b..7761b7554 100644
--- a/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py
+++ b/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py
@@ -4,319 +4,143 @@
 
 import openmdao.api as om
 from openmdao.utils.assert_utils import assert_near_equal
-from openmdao.utils.testing_utils import use_tempdirs
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 import dymos as dm
 from dymos.utils.misc import om_version
 from dymos.examples.aircraft_steady_flight.aircraft_ode import AircraftODE
 
-optimizer = os.environ.get('DYMOS_DEFAULT_OPT', 'SLSQP')
 
 
 @use_tempdirs
+@require_pyoptsparse(optimizer='IPOPT')
 class TestAircraftCruise(unittest.TestCase):
 
-    def test_cruise_results_gl(self):
-        p = om.Problem()
-        if optimizer == 'SNOPT':
-            p.driver = om.pyOptSparseDriver()
-            p.driver.options['optimizer'] = optimizer
-            p.driver.declare_coloring()
-            p.driver.opt_settings['Major iterations limit'] = 100
-            p.driver.opt_settings['Major step limit'] = 0.05
-            p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
-            p.driver.opt_settings['Major optimality tolerance'] = 1.0E-6
-            p.driver.opt_settings["Linesearch tolerance"] = 0.10
-            p.driver.opt_settings['iSumm'] = 6
-            p.driver.opt_settings['Verify level'] = 3
-        else:
-            p.driver = om.ScipyOptimizeDriver()
-            p.driver.declare_coloring()
-
-        # Pass Reference Area from an external source
-        assumptions = p.model.add_subsystem('assumptions', om.IndepVarComp())
-        assumptions.add_output('S', val=427.8, units='m**2')
-        assumptions.add_output('mass_empty', val=1.0, units='kg')
-        assumptions.add_output('mass_payload', val=1.0, units='kg')
-
-        transcription = dm.GaussLobatto(num_segments=1,
-                                        order=13,
-                                        compressed=False)
-        phase = dm.Phase(ode_class=AircraftODE, transcription=transcription)
-        p.model.add_subsystem('phase0', phase)
-
-        phase.set_time_options(fix_initial=True,
-                               duration_bounds=(3600, 3600),
-                               duration_ref=3600)
-
-        phase.add_state('range', units='km', fix_initial=True, fix_final=False, scaler=0.01,
-                        rate_source='range_rate_comp.dXdt:range',
-                        defect_scaler=0.01)
-        phase.add_state('mass_fuel', units='kg', fix_final=True, upper=20000.0, lower=0.0,
-                        rate_source='propulsion.dXdt:mass_fuel',
-                        scaler=1.0E-4, defect_scaler=1.0E-2)
-        phase.add_state('alt',
-                        rate_source='climb_rate',
-                        units='km', fix_initial=True)
-
-        phase.add_control('mach', targets=['tas_comp.mach', 'aero.mach'], units=None, opt=False)
-
-        phase.add_control('climb_rate', targets=['gam_comp.climb_rate'], units='m/s', opt=False)
-
-        phase.add_parameter('S',
-                            targets=['aero.S', 'flight_equilibrium.S', 'propulsion.S'],
-                            units='m**2')
-
-        phase.add_parameter('mass_empty', targets=['mass_comp.mass_empty'], units='kg')
-        phase.add_parameter('mass_payload', targets=['mass_comp.mass_payload'], units='kg')
-
-        phase.add_path_constraint('propulsion.tau', lower=0.01, upper=1.0)
-
-        phase.add_timeseries_output('tas_comp.TAS')
-
-        p.model.connect('assumptions.S', 'phase0.parameters:S')
-        p.model.connect('assumptions.mass_empty', 'phase0.parameters:mass_empty')
-        p.model.connect('assumptions.mass_payload', 'phase0.parameters:mass_payload')
-
-        phase.add_objective('time', loc='final', ref=3600)
-
-        p.model.linear_solver = om.DirectSolver()
-
-        p.setup()
-
-        phase.set_time_val(initial=0.0, duration=3600)
-        phase.set_state_val('range', (0, 1296.4))
-        phase.set_state_val('mass_fuel', (12236.594555, 1))
-        phase.set_state_val('alt', 5.0)
-        phase.set_control_val('mach', 0.8)
-        phase.set_control_val('climb_rate', 0.0)
-
-        p['assumptions.S'] = 427.8
-        p['assumptions.mass_empty'] = 0.15E6
-        p['assumptions.mass_payload'] = 84.02869 * 400
-
-        dm.run_problem(p)
-
-        time = p.get_val('phase0.timeseries.time')
-        tas = p.get_val('phase0.timeseries.TAS', units='km/s')
-        range = p.get_val('phase0.timeseries.range')
-
-        assert_near_equal(range, tas*time, tolerance=1.0E-4)
-
-        exp_out = phase.simulate()
-
-        time = exp_out.get_val('phase0.timeseries.time')
-        tas = exp_out.get_val('phase0.timeseries.TAS', units='km/s')
-        range = exp_out.get_val('phase0.timeseries.range')
-
-        assert_near_equal(range, tas*time, tolerance=1.0E-4)
-
-    def test_cruise_results_radau(self):
-        p = om.Problem(model=om.Group())
-        if optimizer == 'SNOPT':
-            p.driver = om.pyOptSparseDriver()
-            p.driver.options['optimizer'] = optimizer
-            p.driver.declare_coloring()
-            p.driver.opt_settings['Major iterations limit'] = 100
-            p.driver.opt_settings['Major step limit'] = 0.05
-            p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
-            p.driver.opt_settings['Major optimality tolerance'] = 1.0E-6
-            p.driver.opt_settings["Linesearch tolerance"] = 0.10
-            p.driver.opt_settings['iSumm'] = 6
-            p.driver.opt_settings['Verify level'] = 3
-        else:
-            p.driver = om.ScipyOptimizeDriver()
-            p.driver.declare_coloring()
-
-        # Pass Reference Area from an external source
-        assumptions = p.model.add_subsystem('assumptions', om.IndepVarComp())
-        assumptions.add_output('S', val=427.8, units='m**2')
-        assumptions.add_output('mass_empty', val=1.0, units='kg')
-        assumptions.add_output('mass_payload', val=1.0, units='kg')
-
-        transcription = dm.Radau(num_segments=1, order=13, compressed=False)
-        phase = dm.Phase(ode_class=AircraftODE, transcription=transcription)
-        p.model.add_subsystem('phase0', phase)
-
-        phase.set_time_options(fix_initial=True,
-                               duration_bounds=(3600, 3600),
-                               duration_ref=3600)
-
-        phase.add_state('range', units='km', fix_initial=True, fix_final=False, scaler=0.01,
-                        rate_source='range_rate_comp.dXdt:range',
-                        defect_scaler=0.01)
-        phase.add_state('mass_fuel', units='kg', fix_final=True, upper=20000.0, lower=0.0,
-                        rate_source='propulsion.dXdt:mass_fuel',
-                        scaler=1.0E-4, defect_scaler=1.0E-2)
-        phase.add_state('alt',
-                        rate_source='climb_rate',
-                        units='km', fix_initial=True)
-
-        phase.add_control('mach', targets=['tas_comp.mach', 'aero.mach'], opt=False)
-
-        phase.add_control('climb_rate', targets=['gam_comp.climb_rate'], units='m/s', opt=False)
-
-        phase.add_parameter('S',
-                            targets=['aero.S', 'flight_equilibrium.S', 'propulsion.S'],
-                            units='m**2')
-
-        phase.add_parameter('mass_empty', targets=['mass_comp.mass_empty'], units='kg')
-        phase.add_parameter('mass_payload', targets=['mass_comp.mass_payload'], units='kg')
-
-        phase.add_path_constraint('propulsion.tau', lower=0.01, upper=1.0)
-
-        phase.add_timeseries_output('tas_comp.TAS')
-
-        p.model.connect('assumptions.S', 'phase0.parameters:S')
-        p.model.connect('assumptions.mass_empty', 'phase0.parameters:mass_empty')
-        p.model.connect('assumptions.mass_payload', 'phase0.parameters:mass_payload')
-
-        phase.add_objective('time', loc='final', ref=3600)
-
-        p.model.linear_solver = om.DirectSolver()
-
-        p.setup()
-
-        phase.set_time_val(initial=0.0, duration=3600)
-        phase.set_state_val('range', (0, 1296.4))
-        phase.set_state_val('mass_fuel', (12236.594555, 0))
-        phase.set_state_val('alt', 5.0)
-        phase.set_control_val('mach', 0.8)
-        phase.set_control_val('climb_rate', 0.0)
-
-        p['assumptions.S'] = 427.8
-        p['assumptions.mass_empty'] = 0.15E6
-        p['assumptions.mass_payload'] = 84.02869 * 400
-
-        dm.run_problem(p)
-
-        time = p.get_val('phase0.timeseries.time')
-        tas = p.get_val('phase0.timeseries.TAS', units='km/s')
-        range = p.get_val('phase0.timeseries.range')
-
-        assert_near_equal(range, tas*time, tolerance=1.0E-4)
-
-        exp_out = phase.simulate()
-
-        time = exp_out.get_val('phase0.timeseries.time')
-        tas = exp_out.get_val('phase0.timeseries.TAS', units='km/s')
-        range = exp_out.get_val('phase0.timeseries.range')
-
-        assert_near_equal(range, tas*time, tolerance=1.0E-4)
-
-    def test_cruise_results_birkhoff(self):
-        p = om.Problem(model=om.Group())
-        optimizer = 'SLSQP'
-        if optimizer == 'SNOPT':
-            p.driver = om.pyOptSparseDriver()
-            p.driver.options['optimizer'] = optimizer
-            p.driver.declare_coloring()
-            p.driver.opt_settings['Major iterations limit'] = 100
-            p.driver.opt_settings['Major step limit'] = 0.05
-            p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
-            p.driver.opt_settings['Major optimality tolerance'] = 1.0E-6
-            p.driver.opt_settings["Linesearch tolerance"] = 0.10
-            p.driver.opt_settings['iSumm'] = 6
-            p.driver.opt_settings['Verify level'] = 3
-        else:
-            p.driver = om.ScipyOptimizeDriver()
-            p.driver.declare_coloring()
-
-        # Pass Reference Area from an external source
-        assumptions = p.model.add_subsystem('assumptions', om.IndepVarComp())
-        assumptions.add_output('S', val=427.8, units='m**2')
-        assumptions.add_output('mass_empty', val=1.0, units='kg')
-        assumptions.add_output('mass_payload', val=1.0, units='kg')
-
-        transcription = dm.Birkhoff(num_nodes=50, grid_type='lgl', solve_segments='forward')
-        phase = dm.Phase(ode_class=AircraftODE, transcription=transcription)
-        traj = dm.Trajectory()
-        traj.add_phase('phase0', phase)
-        p.model.add_subsystem('traj', traj)
-
-        phase.set_time_options(fix_initial=True,
-                               duration_bounds=(3600, 3600),
-                               duration_ref=3600)
-
-        phase.add_state('range', units='km', fix_initial=True, fix_final=False, scaler=0.01,
-                        rate_source='range_rate_comp.dXdt:range',
-                        defect_scaler=0.01)
-        phase.add_state('mass_fuel', units='kg', fix_final=False, upper=20000.0, lower=0.0,
-                        rate_source='propulsion.dXdt:mass_fuel',
-                        scaler=1.0E-4, defect_scaler=1.0E-2)
-        phase.add_state('alt',
-                        rate_source='climb_rate',
-                        units='km', fix_initial=True)
-
-        phase.add_control('mach', targets=['tas_comp.mach', 'aero.mach'], opt=False)
-
-        phase.add_control('climb_rate', targets=['gam_comp.climb_rate'], units='m/s', opt=False)
-
-        phase.add_parameter('S',
-                            targets=['aero.S', 'flight_equilibrium.S', 'propulsion.S'],
-                            units='m**2')
-
-        phase.add_parameter('mass_empty', targets=['mass_comp.mass_empty'], units='kg')
-        phase.add_parameter('mass_payload', targets=['mass_comp.mass_payload'], units='kg')
-
-        phase.add_path_constraint('propulsion.tau', lower=0.01, upper=1.0)
-        phase.add_boundary_constraint('mass_fuel', loc='final', equals=0.)
-        phase.add_boundary_constraint('range_rate_comp.dXdt:range', loc='final', lower=0.)
-        phase.add_boundary_constraint('range_rate_comp.dXdt:range', loc='initial', lower=0.)
-
-        phase.add_timeseries_output('tas_comp.TAS')
-
-        p.model.connect('assumptions.S', 'traj.phase0.parameters:S')
-        p.model.connect('assumptions.mass_empty', 'traj.phase0.parameters:mass_empty')
-        p.model.connect('assumptions.mass_payload', 'traj.phase0.parameters:mass_payload')
-
-        phase.add_objective('time', loc='final', ref=3600)
-
-        p.model.linear_solver = om.DirectSolver()
-
-        p.setup()
-
-        p['traj.phase0.t_initial'] = 0.0
-        p['traj.phase0.t_duration'] = 1.0 * 3600.0
-
-        p['traj.phase0.initial_states:range'] = 0.0
-        p['traj.phase0.initial_states:mass_fuel'] = 10_000.
-        p['traj.phase0.initial_states:alt'] = 10.0
-
-        p['traj.phase0.final_states:mass_fuel'] = 1.0
-        p['traj.phase0.final_states:range'] = 500.0
-        p['traj.phase0.final_states:alt'] = 10.0
-
-        p['traj.phase0.controls:mach'] = 0.8
-        p['traj.phase0.controls:climb_rate'] = 0.0
-
-        p['assumptions.S'] = 427.8
-        p['assumptions.mass_empty'] = 0.15E6
-        p['assumptions.mass_payload'] = 84.02869 * 400
-
-        dm.run_problem(p, simulate=True, make_plots=True)
-
-        time = p.get_val('traj.phase0.timeseries.time')
-        tas = p.get_val('traj.phase0.timeseries.TAS', units='km/s')
-        range = p.get_val('traj.phase0.timeseries.range')
-        mass_fuel = p.get_val('traj.phase0.timeseries.mass_fuel')
-
-        assert_near_equal(range, tas*time, tolerance=1.0E-4)
-        assert_near_equal(mass_fuel[-1, ...], 0.0, tolerance=1.0E-4)
-
-        if om_version()[0] > (3, 34, 2):
-            sim_out_dir = traj.sim_prob.get_outputs_dir()
-        else:
-            sim_out_dir = pathlib.Path.cwd()
-        case = om.CaseReader(sim_out_dir / 'dymos_simulation.db').get_case('final')
-
-        time = case.get_val('traj.phase0.timeseries.time')
-        tas = case.get_val('traj.phase0.timeseries.TAS', units='km/s')
-        range = case.get_val('traj.phase0.timeseries.range')
-        mass_fuel = case.get_val('traj.phase0.timeseries.mass_fuel')
-
-        assert_near_equal(range, tas*time, tolerance=1.0E-4)
-        assert_near_equal(mass_fuel[-1, ...], 0.0, tolerance=1.0E-4)
+    def test_cruise_results(self):
+
+        txs = {'gauss-lobatto': dm.GaussLobatto(num_segments=1,
+                                                order=13,
+                                                compressed=False),
+               'gauss-lobatto-compressed': dm.GaussLobatto(num_segments=2,
+                                                           order=7,
+                                                           compressed=True),
+               'radau': dm.Radau(num_segments=2,
+                                 order=13,
+                                 compressed=False),
+               'radau-compressed': dm.Radau(num_segments=2,
+                                            order=7,
+                                            compressed=True),
+               'birkhoff': dm.Birkhoff(num_nodes=13,
+                                       compressed=False)}
+
+        for tx_name, tx in txs.items():
+            with self.subTest(msg=f'{tx_name=}'):
+                p = om.Problem(model=om.Group())
+                optimizer = 'IPOPT'
+                if optimizer == 'SNOPT':
+                    p.driver = om.pyOptSparseDriver()
+                    p.driver.options['optimizer'] = optimizer
+                    p.driver.opt_settings['Major iterations limit'] = 100
+                    # p.driver.opt_settings['Major step limit'] = 0.05
+                    p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-8
+                    p.driver.opt_settings['Major optimality tolerance'] = 1.0E-8
+                    p.driver.opt_settings["Linesearch tolerance"] = 0.10
+                    p.driver.opt_settings['iSumm'] = 6
+                elif optimizer == 'IPOPT':
+                    p.driver = om.pyOptSparseDriver(optimizer=optimizer)
+                    p.driver.opt_settings['print_level'] = 5
+                else:
+                    p.driver = om.ScipyOptimizeDriver()
+                p.driver.declare_coloring()
+
+                # Pass Reference Area from an external source
+                assumptions = p.model.add_subsystem('assumptions', om.IndepVarComp())
+                assumptions.add_output('S', val=427.8, units='m**2')
+                assumptions.add_output('mass_empty', val=1.0, units='kg')
+                assumptions.add_output('mass_payload', val=1.0, units='kg')
+
+                # transcription = dm.Birkhoff(num_nodes=50)
+                # transcription = dm.Radau(order=3, num_segments=10, compressed=False)
+                phase = dm.Phase(ode_class=AircraftODE, transcription=tx)
+                traj = dm.Trajectory()
+                traj.add_phase('phase0', phase)
+                p.model.add_subsystem('traj', traj)
+
+                phase.set_time_options(fix_initial=True,
+                                    duration_bounds=(3600, 3600),
+                                    duration_ref=3600)
+
+                phase.add_state('range', units='km', fix_initial=True, fix_final=False, scaler=0.01,
+                                rate_source='range_rate_comp.dXdt:range',
+                                defect_scaler=0.01)
+                phase.add_state('mass_fuel', units='kg', fix_final=True, upper=20000.0, lower=0.0,
+                                rate_source='propulsion.dXdt:mass_fuel',
+                                scaler=1.0E-4, defect_scaler=1.0E-2)
+                phase.add_state('alt',
+                                rate_source='climb_rate',
+                                units='km', fix_initial=True)
+
+                phase.add_control('mach', targets=['tas_comp.mach', 'aero.mach'], opt=False)
+
+                phase.add_control('climb_rate', targets=['gam_comp.climb_rate'], units='m/s', opt=False)
+
+                phase.add_parameter('S',
+                                    targets=['aero.S', 'flight_equilibrium.S', 'propulsion.S'],
+                                    units='m**2')
+
+                phase.add_parameter('mass_empty', targets=['mass_comp.mass_empty'], units='kg')
+                phase.add_parameter('mass_payload', targets=['mass_comp.mass_payload'], units='kg')
+
+                phase.add_path_constraint('propulsion.tau', lower=0.01, upper=1.0)
+
+                phase.add_timeseries_output('tas_comp.TAS')
+
+                p.model.connect('assumptions.S', 'traj.phase0.parameters:S')
+                p.model.connect('assumptions.mass_empty', 'traj.phase0.parameters:mass_empty')
+                p.model.connect('assumptions.mass_payload', 'traj.phase0.parameters:mass_payload')
+
+                phase.add_objective('time', loc='final', ref=3600)
+
+                p.model.linear_solver = om.DirectSolver()
+
+                p.setup()
+
+                phase.set_time_val(initial=0, duration=3600.0)
+
+                phase.set_state_val('range', [0.0, 500.0])
+                phase.set_state_val('mass_fuel', [10_000.0, 0.0])
+                phase.set_state_val('alt', [10.0, 10.0])
+
+                phase.set_control_val('mach', 0.8)
+                phase.set_control_val('climb_rate', 0.0)
+
+                p.set_val('assumptions.S', 427.8)
+                p.set_val('assumptions.mass_empty', 0.15E6)
+                p.set_val('assumptions.mass_payload', 84.02869 * 400)
+
+                dm.run_problem(p, simulate=True)
+
+                time = p.get_val('traj.phase0.timeseries.time')
+                tas = p.get_val('traj.phase0.timeseries.TAS', units='km/s')
+                r = p.get_val('traj.phase0.timeseries.range')
+                mass_fuel = p.get_val('traj.phase0.timeseries.mass_fuel')
+
+                assert_near_equal(r, tas*time, tolerance=1.0E-4)
+                assert_near_equal(mass_fuel[-1, ...], 0.0, tolerance=1.0E-4)
+
+                if om_version()[0] > (3, 34, 2):
+                    sim_out_dir = traj.sim_prob.get_outputs_dir()
+                else:
+                    sim_out_dir = pathlib.Path.cwd()
+                case = om.CaseReader(sim_out_dir / 'dymos_simulation.db').get_case('final')
+
+                time = case.get_val('traj.phase0.timeseries.time')
+                tas = case.get_val('traj.phase0.timeseries.TAS', units='km/s')
+                r = case.get_val('traj.phase0.timeseries.range')
+                mass_fuel = case.get_val('traj.phase0.timeseries.mass_fuel')
+
+                assert_near_equal(r, tas*time, tolerance=1.0E-4)
+                assert_near_equal(mass_fuel[-1, ...], 0.0, tolerance=1.0E-4)
 
 
 if __name__ == '__main__':  # pragma: no cover
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 e4ccd25c7..187be4903 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
@@ -257,3 +257,7 @@ def test_balanced_field_length_for_docs(self):
         assert_near_equal(2114.387, sol_r_f_rto, tolerance=0.01)
         assert_near_equal(2114.387, sim_r_f_climb, tolerance=0.01)
         assert_near_equal(2114.387, sim_r_f_rto, tolerance=0.01)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py
index 6424cbfa6..f95fc4c21 100644
--- a/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py
+++ b/dymos/examples/brachistochrone/doc/test_doc_brachistochrone_static_gravity.py
@@ -97,7 +97,7 @@ def test_brachistochrone_static_gravity(self):
                         fix_initial=True, fix_final=False, solve_segments=False)
 
         phase.add_control('theta', targets=['theta'],
-                          continuity=True, rate_continuity=True,
+                          continuity=False, rate_continuity=False,
                           units='deg', lower=0.01, upper=179.9)
 
         phase.add_parameter('g', targets=['g'], static_target=True, opt=False)
diff --git a/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py b/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py
index acbdcf691..a232c9b34 100644
--- a/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py
+++ b/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py
@@ -51,7 +51,7 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
     # can't fix final position if you're solving the segments
 
     phase.add_state('pos', fix_initial=True, fix_final=fix_final,
-                    solve_segments=solve_segments, ref=[1, 1])
+                    solve_segments=solve_segments)
     #
     phase.add_state('v', fix_initial=True, fix_final=False, solve_segments=solve_segments)
     #
@@ -67,7 +67,7 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
     # Minimize time at the end of the phase
     phase.add_objective('time', loc='final', scaler=10)
 
-    p.model.linear_solver = om.DirectSolver()
+    p.model.linear_solver = om.DirectSolver(assemble_jac=False)
     p.setup(check=True, force_alloc_complex=force_alloc_complex)
 
     phase.set_time_val(initial=0.0, duration=1.8016)
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
index caa52d4df..58d05ddce 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
@@ -138,7 +138,7 @@ def test_brachistochrone_control_rate_targets(self):
 
                     with self.subTest(f'{tx_name=} {control_target_method=} {control_type=}'):
 
-                        p = om.Problem(model=om.Group())
+                        p = om.Problem(model=om.Group(), reports=True)
                         p.driver = om.ScipyOptimizeDriver()
                         p.driver.declare_coloring()
 
@@ -170,7 +170,7 @@ def test_brachistochrone_control_rate_targets(self):
 
                         phase.add_control('int_theta', lower=0.0, upper=None, fix_initial=True,
                                           rate_targets=_unspecified if control_target_method == 'implicit' else ['theta'],
-                                          control_type=control_type, order=7)
+                                          control_type=control_type, order=7, continuity=True)
 
                         phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
 
@@ -187,7 +187,7 @@ def test_brachistochrone_control_rate_targets(self):
 
                         phase.set_state_val('x', [0, 10])
                         phase.set_state_val('y', [10, 5])
-                        phase.set_state_val('v', [0, 9.9])
+                        phase.set_state_val('v', [0.001, 9.9])
                         phase.set_control_val('int_theta', [0, 100], units='deg*s')
 
                         p.run_model()
diff --git a/dymos/examples/brachistochrone/test/test_ex_brachistochrone.py b/dymos/examples/brachistochrone/test/test_ex_brachistochrone.py
index bb6f0fb8f..e785677e4 100644
--- a/dymos/examples/brachistochrone/test/test_ex_brachistochrone.py
+++ b/dymos/examples/brachistochrone/test/test_ex_brachistochrone.py
@@ -1,5 +1,6 @@
 import os
 import unittest
+
 from numpy.testing import assert_almost_equal
 
 import dymos.examples.brachistochrone.test.ex_brachistochrone as ex_brachistochrone
@@ -9,7 +10,7 @@
 OPT, OPTIMIZER = set_pyoptsparse_opt('SNOPT', fallback=True)
 
 
-@use_tempdirs
+# @use_tempdirs
 class TestBrachistochroneExample(unittest.TestCase):
 
     @classmethod
diff --git a/dymos/examples/brachistochrone/test/test_ex_brachistochrone_vector_states.py b/dymos/examples/brachistochrone/test/test_ex_brachistochrone_vector_states.py
index f2ef1f443..d65230624 100644
--- a/dymos/examples/brachistochrone/test/test_ex_brachistochrone_vector_states.py
+++ b/dymos/examples/brachistochrone/test/test_ex_brachistochrone_vector_states.py
@@ -6,11 +6,7 @@
 
 import sys
 
-from openmdao.utils.general_utils import set_pyoptsparse_opt, printoptions
-from openmdao.utils.testing_utils import use_tempdirs, set_env_vars_context
-from openmdao.utils.tests.test_hooks import hooks_active
-
-import dymos as dm
+import dymos
 import dymos.examples.brachistochrone.test.ex_brachistochrone_vector_states as ex_brachistochrone_vs
 from dymos.utils.testing_utils import assert_check_partials, _get_reports_dir
 
@@ -62,36 +58,45 @@ def assert_results(self, p):
 
     def assert_partials(self, p):
         with printoptions(linewidth=1024, edgeitems=100):
-            cpd = p.check_partials(method='cs', compact_print=True, out_stream=None)
+            cpd = p.check_partials(method='cs', compact_print=True,
+                                   show_only_incorrect=True)#, out_stream=None)
         assert_check_partials(cpd)
+        # p.check_totals(method='cs', compact_print=False)
 
     def test_ex_brachistochrone_vs_radau_compressed(self):
         ex_brachistochrone_vs.SHOW_PLOTS = True
         p = ex_brachistochrone_vs.brachistochrone_min_time(transcription='radau-ps',
                                                            compressed=True,
                                                            force_alloc_complex=True,
-                                                           run_driver=False)
+                                                           run_driver=True)
         p.run_driver()
         self.assert_partials(p)
 
     def test_ex_brachistochrone_vs_radau_uncompressed(self):
         ex_brachistochrone_vs.SHOW_PLOTS = True
-        p = ex_brachistochrone_vs.brachistochrone_min_time(transcription='radau-ps',
-                                                           compressed=False,
-                                                           force_alloc_complex=True,
-                                                           run_driver=True)
-        self.assert_results(p)
-        self.assert_partials(p)
+        with dymos.options.temporary(include_check_partials=True):
+            p = ex_brachistochrone_vs.brachistochrone_min_time(transcription='radau-ps',
+                                                               num_segments=3,
+                                                               transcription_order=3,
+                                                               compressed=False,
+                                                               force_alloc_complex=True,
+                                                               dynamic_simul_derivs=False,
+                                                               run_driver=False)
+            # self.assert_results(p)
+            import numpy as np
+            with np.printoptions(linewidth=100_000, edgeitems=100_000):
+                self.assert_partials(p)
 
     def test_ex_brachistochrone_vs_gl_compressed(self):
         ex_brachistochrone_vs.SHOW_PLOTS = True
-        p = ex_brachistochrone_vs.brachistochrone_min_time(transcription='gauss-lobatto',
-                                                           compressed=True,
-                                                           force_alloc_complex=True,
-                                                           run_driver=True)
-
-        self.assert_results(p)
-        self.assert_partials(p)
+        with dymos.options.temporary(include_check_partials=True):
+            p = ex_brachistochrone_vs.brachistochrone_min_time(transcription='gauss-lobatto',
+                                                            compressed=True,
+                                                            force_alloc_complex=True,
+                                                            run_driver=True)
+
+            self.assert_results(p)
+            self.assert_partials(p)
 
     def test_ex_brachistochrone_vs_gl_uncompressed(self):
         ex_brachistochrone_vs.SHOW_PLOTS = True
diff --git a/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py b/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py
index c0d67ac31..d480309f6 100644
--- a/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py
+++ b/dymos/examples/finite_burn_orbit_raise/finite_burn_orbit_raise_problem.py
@@ -287,9 +287,6 @@ def two_burn_orbit_raise_problem(transcription='gauss-lobatto', optimizer='SLSQP
             burn2.set_control_val('u1', [0, 0])
 
     p.final_setup()
-    if burn2 in p.model.traj.phases._subsystems_myproc:
-        print(burn2.get_val('t_initial'))
-        print(burn2.get_val('t_duration'))
 
     if run_driver or simulate:
         dm.run_problem(p, run_driver=run_driver, simulate=simulate, restart=restart, make_plots=True,
diff --git a/dymos/phase/analytic_phase.py b/dymos/phase/analytic_phase.py
index 8c67d8d67..0e3badd0c 100644
--- a/dymos/phase/analytic_phase.py
+++ b/dymos/phase/analytic_phase.py
@@ -124,7 +124,7 @@ def add_control(self, name, units=_unspecified, desc=_unspecified, opt=_unspecif
             This option is invalid if opt=False.
         targets : Sequence of str or None
             Targets in the ODE to which this control is connected.
-            In the future, if left _unspecified (the default), the phase control will try to connect to an ODE input
+            If left _unspecified (the default), the phase control will try to connect to an ODE input
             of the same name. Set targets to None to prevent this.
         rate_targets : Sequence of str or None
             The targets in the ODE to which the control rate is connected.
@@ -212,7 +212,7 @@ def set_control_options(self, name, units=_unspecified, desc=_unspecified, opt=_
             This option is invalid if opt=False.
         targets : Sequence of str or None
             Targets in the ODE to which this control is connected.
-            In the future, if left _unspecified (the default), the phase control will try to connect to an ODE input
+            If left _unspecified (the default), the phase control will try to connect to an ODE input
             of the same name. Set targets to None to prevent this.
         rate_targets : Sequence of str or None
             The targets in the ODE to which the control rate is connected.
diff --git a/dymos/phase/options.py b/dymos/phase/options.py
index 2f3306f40..3a9ef7190 100644
--- a/dymos/phase/options.py
+++ b/dymos/phase/options.py
@@ -506,6 +506,9 @@ def __init__(self, read_only=False):
         self.declare(name='t_duration_targets', allow_none=True, default=[],
                      desc='targets in the ODE to which the total duration of the phase is connected')
 
+        self.declare(name='t_final_targets', allow_none=True, default=[],
+                     desc='targets in the ODE to which the final time of the phase is connected')
+
         self.declare(name='t_duration_balance_options', default={},
                      desc='options dictionary for the duration residual')
 
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index b83fa50ee..979716a21 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -254,7 +254,7 @@ def add_state(self, name, units=_unspecified, shape=_unspecified,
         targets : str or Sequence of str
             The path to the targets of the state variable in the ODE system.  If given
             this will override the value given by the @declare_state decorator on the ODE.
-            In the future, if left _unspecified (the default), the phase variable will try to connect to an ODE input
+            If left _unspecified (the default), the phase variable will try to connect to an ODE input
             of the same name. Set targets to None to prevent this.
         val :  ndarray
             The default value of the state at the state discretization nodes of the phase.
@@ -353,7 +353,7 @@ def set_state_options(self, name, units=_unspecified, shape=_unspecified,
         targets : str or Sequence of str
             The path to the targets of the state variable in the ODE system.  If given
             this will override the value given by the @declare_state decorator on the ODE.
-            In the future, if left _unspecified (the default), the phase variable will try to connect to an ODE input
+            If left _unspecified (the default), the phase variable will try to connect to an ODE input
             of the same name. Set targets to None to prevent this.
         val :  ndarray
             The default value of the state at the state discretization nodes of the phase.
@@ -556,7 +556,7 @@ def add_control(self, name, order=_unspecified, units=_unspecified, desc=_unspec
             This option is invalid if opt=False.
         targets : Sequence of str or None
             Targets in the ODE to which this control is connected.
-            In the future, if left _unspecified (the default), the phase control will try to connect to an ODE input
+            If left _unspecified (the default), the phase control will try to connect to an ODE input
             of the same name. Set targets to None to prevent this.
         rate_targets : Sequence of str or None
             The targets in the ODE to which the control rate is connected.
@@ -677,7 +677,7 @@ def set_control_options(self, name, order=_unspecified, units=_unspecified, desc
             This option is invalid if opt=False.
         targets : Sequence of str or None
             Targets in the ODE to which this control is connected.
-            In the future, if left _unspecified (the default), the phase control will try to connect to an ODE input
+            If left _unspecified (the default), the phase control will try to connect to an ODE input
             of the same name. Set targets to None to prevent this.
         rate_targets : Sequence of str or None
             The targets in the ODE to which the control rate is connected.
@@ -1077,7 +1077,7 @@ def add_parameter(self, name, val=_unspecified, units=_unspecified, opt=False,
             The unit-reference value of the parameter for the optimizer.
         targets : Sequence of str or None
             Targets in the ODE to which this parameter is connected.
-            In the future, if left _unspecified (the default), the phase parameter will try to connect to an ODE input
+            If left _unspecified (the default), the phase parameter will try to connect to an ODE input
             of the same name. Set targets to None to prevent this.
         shape : Sequence of int
             The shape of the parameter.
@@ -1147,7 +1147,7 @@ def set_parameter_options(self, name, val=_unspecified, units=_unspecified, opt=
             The unit-reference value of the parameter for the optimizer.
         targets : Sequence of str or None
             Targets in the ODE to which this parameter is connected.
-            In the future, if left _unspecified (the default), the phase parameter will try to connect to an ODE input
+            If left _unspecified (the default), the phase parameter will try to connect to an ODE input
             of the same name. Set targets to None to prevent this.
         shape : Sequence of int
             The shape of the parameter.
diff --git a/dymos/phase/test/test_time_targets.py b/dymos/phase/test/test_time_targets.py
index 48ce3f768..e785dace8 100644
--- a/dymos/phase/test/test_time_targets.py
+++ b/dymos/phase/test/test_time_targets.py
@@ -136,18 +136,11 @@ def _make_problem(self, transcription, num_seg, transcription_order=3, input_ini
 
         p.setup(check=True, force_alloc_complex=True)
 
-        p['phase0.t_initial'] = 1.0
-        p['phase0.t_duration'] = 2.0
-
-        if transcription == 'explicit-shooting':
-            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['phase0.controls:theta'] = phase.interp('theta', [0.01, 100.5])
+        phase.set_time_val(initial=1.0, duration=2.0)
+        phase.set_state_val('x', vals=[0, 10])
+        phase.set_state_val('y', vals=[10, 5])
+        phase.set_state_val('v', vals=[1.0E-6, 9.9])
+        phase.set_control_val('theta', vals=[0.01, 100.5], units='deg')
 
         return p
 
@@ -226,15 +219,11 @@ def test_radau(self):
 
                 time_phase_all = p[f'phase0.timeseries.{time_name}_phase'].ravel()
 
-                assert_near_equal(p['phase0.rhs_all.time_phase'][-1], 1.8016, tolerance=1.0E-3)
+                assert_near_equal(p['phase0.ode_all.time_phase'][-1], 1.8016, tolerance=1.0E-3)
 
-                assert_near_equal(p['phase0.rhs_all.t_initial'], p['phase0.t_initial'])
+                assert_near_equal(p['phase0.ode_all.time_phase'], time_phase_all)
 
-                assert_near_equal(p['phase0.rhs_all.t_duration'], p['phase0.t_duration'])
-
-                assert_near_equal(p['phase0.rhs_all.time_phase'], time_phase_all)
-
-                assert_near_equal(p['phase0.rhs_all.time'], time_all)
+                assert_near_equal(p['phase0.ode_all.time'], time_all)
 
                 exp_out = p.model.phase0.simulate()
 
@@ -360,15 +349,11 @@ def test_radau_targets_are_inputs(self):
 
         time_phase_all = p['phase0.t_phase']
 
-        assert_near_equal(p['phase0.rhs_all.time_phase'][-1], 1.8016, tolerance=1.0E-3)
-
-        assert_near_equal(p['phase0.rhs_all.t_initial'], p['phase0.t_initial'])
-
-        assert_near_equal(p['phase0.rhs_all.t_duration'], p['phase0.t_duration'])
+        assert_near_equal(p['phase0.ode_all.time_phase'][-1], 1.8016, tolerance=1.0E-3)
 
-        assert_near_equal(p['phase0.rhs_all.time_phase'], time_phase_all)
+        assert_near_equal(p['phase0.ode_all.time_phase'], time_phase_all)
 
-        assert_near_equal(p['phase0.rhs_all.time'], time_all)
+        assert_near_equal(p['phase0.ode_all.time'], time_all)
 
         exp_out = p.model._get_subsystem('phase0').simulate()
 
diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index 4ef4c8bef..5d8936477 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -2,6 +2,7 @@
 import warnings
 
 import numpy as np
+from numpy.testing import assert_allclose
 
 from scipy.interpolate import interp1d
 
@@ -197,13 +198,11 @@ def test_timeseries_radau(self, test_smaller_timeseries=False):
         t_sol = p.get_val('phase0.timeseries.time')
         t_sim = exp_out.get_val('phase0.timeseries.time')
 
-        p.model.list_outputs()
-
         for state_name, rate_name in (('x', 'xdot'), ('y', 'ydot'), ('v', 'vdot')):
             rate_sol = p.get_val(f'phase0.timeseries.{rate_name}')
             rate_sim = exp_out.get_val(f'phase0.timeseries.{rate_name}')
             rate_t_sim = interp1d(t_sim.ravel(), rate_sim.ravel())
-            assert_near_equal(rate_t_sim(t_sol), rate_sol, tolerance=1.0E-3)
+            assert_allclose(rate_t_sim(t_sol), rate_sol, rtol=1.0E-3, atol=1.0E-3)
 
     def test_timeseries_radau_smaller_timeseries(self):
         self.test_timeseries_radau(test_smaller_timeseries=True)
diff --git a/dymos/transcriptions/__init__.py b/dymos/transcriptions/__init__.py
index d4e2a1fc1..00b8ff7fb 100644
--- a/dymos/transcriptions/__init__.py
+++ b/dymos/transcriptions/__init__.py
@@ -1,5 +1,6 @@
 from .analytic.analytic import Analytic
 from .explicit_shooting import ExplicitShooting
 from .pseudospectral.gauss_lobatto import GaussLobatto
-from .pseudospectral.radau_pseudospectral import Radau
+from .pseudospectral.radau_pseudospectral import Radau as RadauDeprecated
 from .pseudospectral.birkhoff import Birkhoff
+from .pseudospectral.radau_new import RadauNew as Radau
\ No newline at end of file
diff --git a/dymos/transcriptions/analytic/analytic.py b/dymos/transcriptions/analytic/analytic.py
index d063d05d7..fa5d14fb1 100644
--- a/dymos/transcriptions/analytic/analytic.py
+++ b/dymos/transcriptions/analytic/analytic.py
@@ -85,7 +85,8 @@ def configure_time(self, phase):
                               flat_src_indices=True if dynamic else None)
 
         for name, targets in [('t_initial', options['t_initial_targets']),
-                              ('t_duration', options['t_duration_targets'])]:
+                              ('t_duration', options['t_duration_targets']),
+                              ('t_final', options['t_final_targets'])]:
             for t in targets:
                 shape = ode_inputs[t]['shape']
 
diff --git a/dymos/transcriptions/common/__init__.py b/dymos/transcriptions/common/__init__.py
index 273402be4..318948e9b 100644
--- a/dymos/transcriptions/common/__init__.py
+++ b/dymos/transcriptions/common/__init__.py
@@ -1,5 +1,6 @@
 from .continuity_comp import RadauPSContinuityComp, GaussLobattoContinuityComp
-from .control_group import ControlGroup
+from .control_group import ControlGroup, ControlInterpComp
+from .control_comp import ControlComp
 from .parameter_comp import ParameterComp
 from .time_comp import TimeComp
 from .timeseries_group import TimeseriesOutputGroup
diff --git a/dymos/transcriptions/common/control_comp.py b/dymos/transcriptions/common/control_comp.py
new file mode 100644
index 000000000..054c0328a
--- /dev/null
+++ b/dymos/transcriptions/common/control_comp.py
@@ -0,0 +1,574 @@
+import numpy as np
+import scipy.sparse as sp
+from scipy.linalg import block_diag
+
+import openmdao.api as om
+
+from ..grid_data import GridData
+from ...utils.misc import get_rate_units
+from ...utils.constants import INF_BOUND
+from ...utils.lgl import lgl
+from ...utils.lagrange import lagrange_matrices
+from ..._options import options as dymos_options
+
+
+class ControlComp(om.ExplicitComponent):
+    r"""
+    Class definition for the ControlComp.
+
+    Compute the approximated control values and rates given the values of a control at all nodes,
+    given values at the control discretization nodes.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+
+    Notes
+    -----
+    .. math::
+
+        u = \left[ L \right] u_d
+
+        \dot{u} = \frac{d\tau_s}{dt} \left[ D \right] u_d
+
+        \ddot{u} = \left( \frac{d\tau_s}{dt} \right)^2 \left[ D_2 \right] u_d
+
+    where
+    :math:`u_d` are the values of the control at the control discretization nodes,
+    :math:`u` are the values of the control at all nodes,
+    :math:`\dot{u}` are the time-derivatives of the control at all nodes,
+    :math:`\ddot{u}` are the second time-derivatives of the control at all nodes,
+    :math:`L` is the Lagrange interpolation matrix,
+    :math:`D` is the Lagrange differentiation matrix,
+    and :math:`\frac{d\tau_s}{dt}` is the ratio of segment duration in segment tau space
+    [-1 1] to segment duration in time.
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._no_check_partials = not dymos_options['include_check_partials']
+
+    def initialize(self):
+        """
+        Declare component options.
+        """
+        self.options.declare(
+            'control_options', types=dict,
+            desc='Dictionary of options for the dynamic controls')
+        self.options.declare(
+            'time_units', default=None, allow_none=True, types=str,
+            desc='Units of time')
+        self.options.declare('compressed', types=bool, default=False,
+                             desc='Specifies whether the controls are giving in a compressed '
+                             'transcription (values are specified once at segment boundaries), '
+                             'or not (each end of the boundary can have a unique value). Controls '
+                             'in the Radau transcription are never compressed.')
+        self.options.declare('grid_data', types=GridData, desc='Container object for grid info for the control inputs.')
+        self.options.declare('output_grid_data', types=GridData, allow_none=True, default=None,
+                             desc='GridData object for the output grid. If None, use the same grid_data as the inputs.')
+
+        # Save the names of the dynamic controls/parameters
+        self._input_names = {}
+        self._output_val_names = {}
+        self._output_rate_names = {}
+        self._output_rate2_names = {}
+        self._output_cnty_defect_names = {}
+        self._output_rate_cnty_defect_names = {}
+        self._output_rate2_cnty_defect_names = {}
+        self._matrices = {}
+
+    def setup(self):
+        """
+        Perform setup procedure for the Control interpolation component.
+        """
+        gd = self.options['grid_data']
+        ogd = self.options['output_grid_data'] or gd
+        if self.options['output_grid_data'] is None:
+            self.options['output_grid_data'] = ogd
+
+        if not gd.is_aligned_with(ogd):
+            raise RuntimeError(f'{self.pathname}: The input grid and the output grid must have the same number of '
+                               f'segments and segment spacing, but the input grid segment ends are '
+                               f'\n{gd.segment_ends}\n and the output grid segment ends are \n'
+                               f'{ogd.segment_ends}.')
+
+    def _configure_controls(self):
+        gd = self.options['grid_data']
+        ogd = self.options['output_grid_data'] or gd
+        num_seg = gd.num_segments
+        segment_end_idxs = gd.subset_node_indices['segment_ends']
+        eval_nodes = ogd.node_ptau
+        control_options = self.options['control_options']
+        num_output_nodes = ogd.num_nodes
+        time_units = self.options['time_units']
+
+        for name, options in control_options.items():
+            if 'control_type' not in options:
+                options['control_type'] = 'full'
+            if options['control_type'] == 'polynomial':
+                disc_nodes, _ = lgl(options['order'] + 1)
+                num_control_input_nodes = len(disc_nodes)
+                shape = options['shape']
+                size = np.prod(shape)
+                units = options['units']
+                rate_units = get_rate_units(units, self.options['time_units'], deriv=1)
+                rate2_units = get_rate_units(units, self.options['time_units'], deriv=2)
+
+                input_shape = (num_control_input_nodes,) + shape
+                output_shape = (num_output_nodes,) + shape
+
+                L_de, D_de = lagrange_matrices(disc_nodes, eval_nodes)
+                _, D_dd = lagrange_matrices(disc_nodes, disc_nodes)
+                D2_de = np.dot(D_de, D_dd)
+
+                self._matrices[name] = L_de, D_de, D2_de
+
+                self._input_names[name] = f'controls:{name}'
+                self._output_val_names[name] = f'control_values:{name}'
+                self._output_rate_names[name] = f'control_rates:{name}_rate'
+                self._output_rate2_names[name] = f'control_rates:{name}_rate2'
+
+                self.add_input(self._input_names[name], val=np.ones(input_shape), units=units)
+                self.add_output(self._output_val_names[name], shape=output_shape, units=units)
+                self.add_output(self._output_rate_names[name], shape=output_shape, units=rate_units)
+                self.add_output(self._output_rate2_names[name], shape=output_shape, units=rate2_units)
+
+                self.val_jacs[name] = np.zeros((num_output_nodes, size, num_control_input_nodes, size))
+                self.rate_jacs[name] = np.zeros((num_output_nodes, size, num_control_input_nodes, size))
+                self.rate2_jacs[name] = np.zeros((num_output_nodes, size, num_control_input_nodes, size))
+
+                for i in range(size):
+                    self.val_jacs[name][:, i, :, i] = L_de
+                    self.rate_jacs[name][:, i, :, i] = D_de
+                    self.rate2_jacs[name][:, i, :, i] = D2_de
+
+                self.val_jacs[name] = self.val_jacs[name].reshape((num_output_nodes * size,
+                                                                  num_control_input_nodes * size),
+                                                                  order='C')
+                self.rate_jacs[name] = self.rate_jacs[name].reshape((num_output_nodes * size,
+                                                                    num_control_input_nodes * size),
+                                                                    order='C')
+                self.rate2_jacs[name] = self.rate2_jacs[name].reshape((num_output_nodes * size,
+                                                                      num_control_input_nodes * size),
+                                                                      order='C')
+                self.val_jac_rows[name], self.val_jac_cols[name] = \
+                    np.where(self.val_jacs[name] != 0)
+                self.rate_jac_rows[name], self.rate_jac_cols[name] = \
+                    np.where(self.rate_jacs[name] != 0)
+                self.rate2_jac_rows[name], self.rate2_jac_cols[name] = \
+                    np.where(self.rate2_jacs[name] != 0)
+
+                self.sizes[name] = size
+
+                rs, cs = self.val_jac_rows[name], self.val_jac_cols[name]
+                self.declare_partials(of=self._output_val_names[name],
+                                      wrt=self._input_names[name],
+                                      rows=rs, cols=cs, val=self.val_jacs[name][rs, cs])
+
+                rs = np.concatenate([np.arange(0, num_output_nodes * size, size, dtype=int) + i
+                                    for i in range(size)])
+
+                self.declare_partials(of=self._output_rate_names[name],
+                                      wrt='t_duration', rows=rs, cols=np.zeros_like(rs))
+
+                self.declare_partials(of=self._output_rate_names[name],
+                                      wrt=self._input_names[name],
+                                      rows=self.rate_jac_rows[name], cols=self.rate_jac_cols[name])
+
+                self.declare_partials(of=self._output_rate2_names[name],
+                                      wrt='t_duration', rows=rs, cols=np.zeros_like(rs))
+
+                self.declare_partials(of=self._output_rate2_names[name],
+                                      wrt=self._input_names[name],
+                                      rows=self.rate2_jac_rows[name], cols=self.rate2_jac_cols[name])
+
+            else:
+                self._input_names[name] = f'controls:{name}'
+                self._output_val_names[name] = f'control_values:{name}'
+                self._output_rate_names[name] = f'control_rates:{name}_rate'
+                self._output_rate2_names[name] = f'control_rates:{name}_rate2'
+                self._output_cnty_defect_names[name] = f'cnty_defect_controls:{name}'
+                self._output_rate_cnty_defect_names[name] = f'cnty_defect_control_rates:{name}_rate'
+                self._output_rate2_cnty_defect_names[name] = f'cnty_defect_control_rates:{name}_rate2'
+
+                shape = options['shape']
+                num_control_input_nodes = gd.subset_num_nodes['control_input']
+                ncinps = gd.subset_num_nodes_per_segment['control_input']
+                input_shape = (num_control_input_nodes,) + shape
+                output_shape = (num_output_nodes,) + shape
+                cnty_output_shape = (num_seg - 1,) + shape
+
+                units = options['units']
+                rate_units = get_rate_units(units, time_units)
+                rate2_units = get_rate_units(units, time_units, deriv=2)
+
+                self.add_input(self._input_names[name], val=np.ones(input_shape), units=units)
+
+                self.add_output(self._output_val_names[name], shape=output_shape, units=units)
+
+                self.add_output(self._output_rate_names[name], shape=output_shape, units=rate_units)
+
+                self.add_output(self._output_rate2_names[name], shape=output_shape,
+                                units=rate2_units)
+
+                if gd.num_segments > 1 and not self.options['compressed']:
+                    if options['continuity']:
+                        self.add_output(self._output_cnty_defect_names[name],
+                                        shape=cnty_output_shape, units=units)
+                        self.add_constraint(name=self._output_cnty_defect_names[name],
+                                            scaler=options['continuity_scaler'],
+                                            ref=options['continuity_ref'],
+                                            equals=0.0, linear=True)
+
+                    if options['rate_continuity']:
+                        self.add_output(self._output_rate_cnty_defect_names[name],
+                                        shape=cnty_output_shape, units=rate_units)
+                        self.add_constraint(name=self._output_rate_cnty_defect_names[name],
+                                            scaler=options['rate_continuity_scaler'],
+                                            ref=options['rate_continuity_ref'],
+                                            equals=0.0, linear=True)
+
+                    if options['rate2_continuity']:
+                        self.add_output(self._output_rate2_cnty_defect_names[name],
+                                        shape=cnty_output_shape, units=rate2_units)
+                        self.add_constraint(name=self._output_rate2_cnty_defect_names[name],
+                                            scaler=options['rate2_continuity_scaler'],
+                                            ref=options['rate2_continuity_ref'],
+                                            equals=0.0, linear=True)
+
+                size = np.prod(shape)
+                self.sizes[name] = size
+                sp_eye = sp.eye(size, format='csr')
+
+                # The partial of interpolated value wrt the control input values is linear
+                # and can be computed as the kronecker product of the interpolation matrix (L)
+                # and eye(size).
+                J_val = sp.kron(self.L, sp_eye, format='csr')
+                rs, cs, data = sp.find(J_val)
+                self.declare_partials(of=self._output_val_names[name],
+                                      wrt=self._input_names[name],
+                                      rows=rs, cols=cs, val=data)
+
+                # The jacobian for continuity defects is similar to the jacobian of the inteprpolated values
+                if options['continuity'] and gd.num_segments > 1 and not self.options['compressed']:
+                    J_cnty_def = -J_val[segment_end_idxs, ...][1:-1:2, ...]
+                    J_cnty_def = J_cnty_def.tolil(copy=True)
+                    one_col_idxs = np.cumsum(np.asarray(ncinps[:-1]))
+                    J_cnty_def[np.arange(num_seg -1, dtype=int), one_col_idxs] = 1.0
+
+                    rs, cs, data = sp.find(J_cnty_def)
+
+                    self.declare_partials(of=self._output_cnty_defect_names[name],
+                                        wrt=self._input_names[name],
+                                        rows=rs, cols=cs, val=data)
+
+                # The partials of the output rate and second derivative wrt dt_dstau
+                rs = np.arange(num_output_nodes * size, dtype=int)
+                cs = np.repeat(np.arange(num_output_nodes, dtype=int), size)
+
+                self.declare_partials(of=self._output_rate_names[name],
+                                      wrt='dt_dstau',
+                                      rows=rs, cols=cs)
+
+                self.declare_partials(of=self._output_rate2_names[name],
+                                      wrt='dt_dstau',
+                                      rows=rs, cols=cs)
+
+                # The partials of the rates and second derivatives are nonlinear but the sparsity
+                # pattern is obtained from the kronecker product of the 1st and 2nd differentiation
+                # matrices (D and D2) and eye(size).
+                self.rate_jacs[name] = sp.kron(self.D, sp_eye, format='csr')
+                rs, cs = self.rate_jacs[name].nonzero()
+
+                self.declare_partials(of=self._output_rate_names[name],
+                                      wrt=self._input_names[name],
+                                      rows=rs, cols=cs)
+
+                # Similar with the continuity rate defects.
+                if options['rate_continuity'] and gd.num_segments > 1 and not self.options['compressed']:
+                    J_rate_cnty_def = sp.kron(self.D, sp_eye, format='lil')[segment_end_idxs, ...]
+                    J_rate_cnty_def = J_rate_cnty_def[1:-1]
+                    J_rate_cnty_def[::2] *= -1
+                    J_rate_cnty_def /= gd.node_dptau_dstau[segment_end_idxs, np.newaxis][1:-1]
+                    # This gives twice as many rows as we want. Each odd-index row needs to
+                    # have its elements moved to the row before it.
+                    rs, cs, data = sp.find(J_rate_cnty_def)
+                    rs = rs // 2
+
+                    self.declare_partials(of=self._output_rate_cnty_defect_names[name],
+                                        wrt=self._input_names[name],
+                                        rows=rs, cols=cs, val=data)
+
+                self.rate2_jacs[name] = sp.kron(self.D2, sp_eye, format='csr')
+
+                rs, cs = self.rate2_jacs[name].nonzero()
+
+                self.declare_partials(of=self._output_rate2_names[name],
+                                      wrt=self._input_names[name],
+                                      rows=rs, cols=cs)
+
+                # Similar with the continuity rate defects.
+                if options['rate2_continuity'] and gd.num_segments > 1 and not self.options['compressed']:
+                    J_rate2_cnty_def = sp.kron(self.D2, sp_eye, format='lil')[segment_end_idxs, ...]
+                    J_rate2_cnty_def = J_rate2_cnty_def[1:-1]
+                    J_rate2_cnty_def[::2] *= -1
+                    J_rate2_cnty_def /= gd.node_dptau_dstau[segment_end_idxs, np.newaxis][1:-1]
+                    # This gives twice as many rows as we want. Each odd-index row needs to
+                    # have its elements moved to the row before it.
+                    rs, cs, data = sp.find(J_rate2_cnty_def)
+                    rs = rs // 2
+
+                    self.declare_partials(of=self._output_rate2_cnty_defect_names[name],
+                                          wrt=self._input_names[name],
+                                          rows=rs, cols=cs, val=data)
+
+            if options['opt']:
+                desvar_indices = np.arange(num_control_input_nodes, dtype=int)
+                if options['fix_initial']:
+                    desvar_indices = desvar_indices[1:]
+                if options['fix_final']:
+                    desvar_indices = desvar_indices[:-1]
+                if not(options['fix_initial'] or options['fix_final']):
+                    desvar_indices = None
+
+                lb = -INF_BOUND if options['lower'] is None else options['lower']
+                ub = INF_BOUND if options['upper'] is None else options['upper']
+
+                self.add_design_var(f'controls:{name}',
+                                    lower=lb,
+                                    upper=ub,
+                                    ref=options['ref'],
+                                    ref0=options['ref0'],
+                                    adder=options['adder'],
+                                    scaler=options['scaler'],
+                                    indices=desvar_indices,
+                                    flat_indices=True)
+
+
+    def configure_io(self):
+        """
+        I/O creation is delayed until configure so we can determine shape and units for the states.
+        """
+        time_units = self.options['time_units']
+        gd = self.options['grid_data']
+        ogd = self.options['output_grid_data'] or gd
+        output_num_nodes = ogd.subset_num_nodes['all']
+        num_seg = gd.num_segments
+
+        self.add_input('dt_dstau', shape=output_num_nodes, units=time_units)
+
+        self.val_jacs = {}
+        self.rate_jacs = {}
+        self.rate2_jacs = {}
+        self.val_jac_rows = {}
+        self.val_jac_cols = {}
+        self.rate_jac_rows = {}
+        self.rate_jac_cols = {}
+        self.rate2_jac_rows = {}
+        self.rate2_jac_cols = {}
+        self.sizes = {}
+
+        self.add_input('t_duration', val=1.0, units=self.options['time_units'],
+                       desc='duration of the phase to which this interpolated control group '
+                            'belongs')
+
+        num_disc_nodes = gd.subset_num_nodes['control_disc']
+        num_input_nodes = gd.subset_num_nodes['control_input']
+
+        # Find the indexing matrix that, multiplied by the values at the input nodes,
+        # gives the values at the discretization nodes
+        L_id = np.zeros((num_disc_nodes, num_input_nodes), dtype=float)
+        L_id[np.arange(num_disc_nodes, dtype=int),
+             gd.input_maps['dynamic_control_input_to_disc']] = 1.0
+        L_id = sp.csr_matrix(L_id)
+
+        # Matrices L_da and D_da interpolate values and rates (respectively) at all nodes from
+        # values specified at control discretization nodes.
+        # If the output grid is different than the input grid, then we have to build these these matrices ourselves.
+        if ogd is gd:
+            L_da, D_da = gd.phase_lagrange_matrices('control_disc', 'all', sparse=True)
+        else:
+            L_blocks = []
+            D_blocks = []
+
+            for iseg in range(num_seg):
+                i1, i2 = gd.subset_segment_indices['control_disc'][iseg, :]
+                indices = gd.subset_node_indices['control_disc'][i1:i2]
+                nodes_given = gd.node_stau[indices]
+
+                i1, i2 = ogd.subset_segment_indices['all'][iseg, :]
+                indices = ogd.subset_node_indices['all'][i1:i2]
+                nodes_eval = ogd.node_stau[indices]
+
+                L_block, D_block = lagrange_matrices(nodes_given, nodes_eval)
+
+                L_blocks.append(L_block)
+                D_blocks.append(D_block)
+
+            L_da = sp.csr_matrix(block_diag(*L_blocks))
+            D_da = sp.csr_matrix(block_diag(*D_blocks))
+
+        self.L = L_da.dot(L_id)
+        self.D = D_da.dot(L_id)
+
+        # Matrix D_dd interpolates rates at discretization nodes from values given at control
+        # discretization nodes.
+        _, D_dd = gd.phase_lagrange_matrices('control_disc', 'control_disc', sparse=True)
+
+        # Matrix D2 provides second derivatives at all nodes given values at input nodes.
+        self.D2 = D_da.dot(D_dd.dot(L_id))
+
+        self._configure_controls()
+
+    def compute(self, inputs, outputs):
+        """
+        Compute interpolated control values and rates.
+
+        Parameters
+        ----------
+        inputs : `Vector`
+            `Vector` containing inputs.
+        outputs : `Vector`
+            `Vector` containing outputs.
+        """
+        gd = self.options['grid_data']
+        ogd = self.options['output_grid_data'] or gd
+        control_options = self.options['control_options']
+        num_control_input_nodes = self.options['grid_data'].subset_num_nodes['control_input']
+        num_output_nodes = ogd.num_nodes
+        seg_end_idxs = gd.subset_node_indices['segment_ends']
+        dt_dptau = 0.5 * inputs['t_duration']
+        dptau_dstau = gd.node_dptau_dstau
+
+        for name, options in control_options.items():
+            if control_options[name]['control_type'] == 'polynomial':
+                L_de, D_de, D2_de = self._matrices[name]
+
+                u = inputs[self._input_names[name]]
+
+                a = np.tensordot(D_de, u, axes=(1, 0)).T
+                b = np.tensordot(D2_de, u, axes=(1, 0)).T
+
+                # divide each "row" by dt_dptau or dt_dptau**2
+                outputs[self._output_val_names[name]] = np.tensordot(L_de, u, axes=(1, 0))
+                outputs[self._output_rate_names[name]] = (a / dt_dptau).T
+                outputs[self._output_rate2_names[name]] = (b / dt_dptau ** 2).T
+            else:
+                size = np.prod(options['shape'])
+
+                u_flat = np.reshape(inputs[self._input_names[name]],
+                                    (num_control_input_nodes, size))
+
+                a = self.D.dot(u_flat)  # du/dstau
+                b = self.D2.dot(u_flat) # d2u/dstau2
+
+                val = np.reshape(self.L.dot(u_flat), (num_output_nodes,) + options['shape'])
+
+                rate = a / inputs['dt_dstau'][:, np.newaxis]
+                rate = np.reshape(rate, (num_output_nodes,) + options['shape'])
+
+                rate2 = b / inputs['dt_dstau'][:, np.newaxis] ** 2
+                rate2 = np.reshape(rate2, (num_output_nodes,) + options['shape'])
+
+                outputs[self._output_val_names[name]] = val
+                outputs[self._output_rate_names[name]] = rate
+                outputs[self._output_rate2_names[name]] = rate2
+
+                if options['continuity'] and gd.num_segments > 1 and not self.options['compressed']:
+                    # output_name = self._output_val_names[name]
+                    cnty_output_name = self._output_cnty_defect_names[name]
+                    segment_end_values = val[seg_end_idxs, ...]
+                    end_vals = segment_end_values[1:-1:2, ...]
+                    start_vals = segment_end_values[2:-1:2, ...]
+                    outputs[cnty_output_name] = start_vals - end_vals
+
+                if options['rate_continuity'] and gd.num_segments > 1 and not self.options['compressed']:
+                    # output_name = self._output_rate_names[name]
+                    rate_in_ptau = a / dptau_dstau[:, np.newaxis]
+                    rate_in_ptau = np.reshape(rate_in_ptau, (num_output_nodes,) + options['shape'])
+
+                    cnty_output_name = self._output_rate_cnty_defect_names[name]
+                    segment_end_values = rate_in_ptau[seg_end_idxs, ...]
+                    end_vals = segment_end_values[1:-1:2, ...]
+                    start_vals = segment_end_values[2:-1:2, ...]
+                    outputs[cnty_output_name] = (start_vals - end_vals)
+
+                if options['rate2_continuity'] and gd.num_segments > 1 and not self.options['compressed']:
+                    # output_name  = self._output_rate2_names[name]
+                    cnty_output_name = self._output_rate2_cnty_defect_names[name]
+                    segment_end_values = rate2[seg_end_idxs, ...]
+                    end_vals = segment_end_values[1:-1:2, ...]
+                    start_vals = segment_end_values[2:-1:2, ...]
+                    outputs[cnty_output_name] = (start_vals - end_vals) #* dt_dptau ** 2
+
+    def compute_partials(self, inputs, partials):
+        """
+        Compute sub-jacobian parts. The model is assumed to be in an unscaled state.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        partials : Jacobian
+            Subjac components written to partials[output_name, input_name].
+        """
+        control_options = self.options['control_options']
+
+        dstau_dt = np.reciprocal(inputs['dt_dstau'])
+        dstau_dt2 = (dstau_dt ** 2)[:, np.newaxis]
+        dstau_dt3 = (dstau_dt ** 3)[:, np.newaxis]
+
+        t_duration = inputs['t_duration']
+
+        for name in control_options:
+            if control_options[name]['control_type'] == 'polynomial':
+                options = control_options[name]
+                control_name = self._input_names[name]
+                num_input_nodes = options['order'] + 1
+                L_de, D_de, D2_de = self._matrices[name]
+                nn = self.options['output_grid_data'].num_nodes
+
+                size = self.sizes[name]
+                rate_name = self._output_rate_names[name]
+                rate2_name = self._output_rate2_names[name]
+
+                # Unroll matrix-shaped controls into an array at each node
+                u_d = np.reshape(inputs[control_name], (num_input_nodes, size))
+
+                t_duration_tile = np.tile(t_duration, size * nn)
+
+                partials[rate_name, 't_duration'] = \
+                    0.5 * (-np.dot(D_de, u_d).ravel(order='F') / (0.5 * t_duration_tile) ** 2)
+
+                partials[rate2_name, 't_duration'] = \
+                    -1.0 * (np.dot(D2_de, u_d).ravel(order='F') / (0.5 * t_duration_tile) ** 3)
+
+                t_duration_x_size = np.repeat(t_duration, size * nn)[:, np.newaxis]
+
+                r_nz, c_nz = self.rate_jac_rows[name], self.rate_jac_cols[name]
+                partials[rate_name, control_name] = \
+                    (self.rate_jacs[name] / (0.5 * t_duration_x_size))[r_nz, c_nz]
+
+                r_nz, c_nz = self.rate2_jac_rows[name], self.rate2_jac_cols[name]
+                partials[rate2_name, control_name] = \
+                    (self.rate2_jacs[name] / (0.5 * t_duration_x_size) ** 2)[r_nz, c_nz]
+            else:
+                control_name = self._input_names[name]
+
+                size = self.sizes[name]
+                rate_name = self._output_rate_names[name]
+                rate2_name = self._output_rate2_names[name]
+                num_input_nodes = self.options['grid_data'].subset_num_nodes['control_input']
+
+                # Unroll shaped controls into an array at each node
+                u_flat = np.reshape(inputs[control_name], (num_input_nodes, size))
+
+                partials[rate_name, 'dt_dstau'] = (-self.D.dot(u_flat) * dstau_dt2).ravel()
+                partials[rate2_name, 'dt_dstau'] = (-2.0 * self.D2.dot(u_flat) * dstau_dt3).ravel()
+
+                dstau_dt_x_size = np.repeat(dstau_dt, size)[:, np.newaxis]
+                dstau_dt2_x_size = np.repeat(dstau_dt2, size)[:, np.newaxis]
+
+                partials[rate_name, control_name] = self.rate_jacs[name].multiply(dstau_dt_x_size).data
+
+                partials[rate2_name, control_name] = self.rate2_jacs[name].multiply(dstau_dt2_x_size).data
\ No newline at end of file
diff --git a/dymos/transcriptions/common/test/test_time_comp.py b/dymos/transcriptions/common/test/test_time_comp.py
index 79e8a7b09..d28d70470 100644
--- a/dymos/transcriptions/common/test/test_time_comp.py
+++ b/dymos/transcriptions/common/test/test_time_comp.py
@@ -1,10 +1,10 @@
 import unittest
 
 import numpy as np
-from numpy.testing import assert_almost_equal
 
 import openmdao.api as om
 from openmdao.utils.testing_utils import use_tempdirs
+from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials
 
 from dymos.transcriptions.grid_data import GridData
 from dymos.transcriptions.common import TimeComp
@@ -35,24 +35,19 @@ def test_results_gauss_lobatto(self):
 
         p = om.Problem(model=om.Group())
 
-        ivc = p.model.add_subsystem(name='ivc', subsys=om.IndepVarComp(), promotes_outputs=['*'])
-
-        ivc.add_output('t_initial', val=0.0, units='s')
-        ivc.add_output('t_duration', val=100.0, units='s')
-
         p.model.add_subsystem('time',
                               TimeComp(num_nodes=gd.num_nodes, node_ptau=gd.node_ptau,
                                        node_dptau_dstau=gd.node_dptau_dstau, units='s'))
 
-        p.model.connect('t_initial', 'time.t_initial')
-        p.model.connect('t_duration', 'time.t_duration')
+        p.setup(check=False, force_alloc_complex=True)
 
-        p.setup(check=True, force_alloc_complex=True)
+        p.set_val('time.t_initial', 0.0, units='s')
+        p.set_val('time.t_duration', 100.0, units='s')
 
         p.run_model()
 
         # Affine transform the nodes [0, 100]
-        dt_dptau = (p['t_duration'] - p['t_initial']) / 2.0
+        dt_dptau = (p['time.t_duration'] - p['time.t_initial']) / 2.0
 
         nodes = []
 
@@ -78,8 +73,12 @@ def test_results_gauss_lobatto(self):
 
             nodes.extend(nodes_time.tolist())
 
-        assert_almost_equal(p['time.t'], nodes)
-        assert_almost_equal(p['time.dt_dstau'], dt_dstau_per_node)
+        assert_near_equal(p['time.t'], nodes)
+        assert_near_equal(p['time.dt_dstau'], dt_dstau_per_node)
+        assert_near_equal(p.get_val('time.t_final'), 100.0)
+
+        cpd = p.check_partials(out_stream=None)
+        assert_check_partials(cpd)
 
     def test_results_radau(self):
 
@@ -90,24 +89,19 @@ def test_results_radau(self):
 
         p = om.Problem(model=om.Group())
 
-        ivc = p.model.add_subsystem(name='ivc', subsys=om.IndepVarComp(), promotes_outputs=['*'])
-
-        ivc.add_output('t_initial', val=0.0, units='s')
-        ivc.add_output('t_duration', val=100.0, units='s')
-
         p.model.add_subsystem('time',
                               TimeComp(num_nodes=gd.num_nodes, node_ptau=gd.node_ptau,
                                        node_dptau_dstau=gd.node_dptau_dstau, units='s'))
 
-        p.model.connect('t_initial', 'time.t_initial')
-        p.model.connect('t_duration', 'time.t_duration')
+        p.setup(check=False, force_alloc_complex=True)
 
-        p.setup(check=True, force_alloc_complex=True)
+        p.set_val('time.t_initial', 0.0, units='s')
+        p.set_val('time.t_duration', 100.0, units='s')
 
         p.run_model()
 
         # Affine transform the nodes [0, 100]
-        dt_dptau = (p['t_duration'] - p['t_initial']) / 2.0
+        dt_dptau = (p['time.t_duration'] - p['time.t_initial']) / 2.0
 
         nodes = []
 
@@ -133,8 +127,12 @@ def test_results_radau(self):
 
             nodes.extend(nodes_time.tolist())
 
-        assert_almost_equal(p['time.t'], nodes, decimal=4)
-        assert_almost_equal(p['time.dt_dstau'], dt_dstau_per_node)
+        assert_near_equal(p['time.t'], nodes, tolerance=1.0E-4)
+        assert_near_equal(p['time.dt_dstau'], dt_dstau_per_node)
+        assert_near_equal(p.get_val('time.t_final'), 100.0)
+
+        cpd = p.check_partials(out_stream=None)
+        assert_check_partials(cpd)
 
 
 if __name__ == "__main__":
diff --git a/dymos/transcriptions/common/time_comp.py b/dymos/transcriptions/common/time_comp.py
index 00bf4a2a9..bfc76fade 100644
--- a/dymos/transcriptions/common/time_comp.py
+++ b/dymos/transcriptions/common/time_comp.py
@@ -57,6 +57,9 @@ def configure_io(self):
         self.add_input('t_duration', val=self.options['duration_val'], units=time_units)
         self.add_output('t', units=time_units, val=np.ones(num_nodes))
         self.add_output('t_phase', units=time_units, val=np.ones(num_nodes))
+        self.add_output('t_final',
+                        val=self.options['initial_val'] + self.options['duration_val'],
+                        units=time_units)
         self.add_output('dt_dstau', units=time_units, val=np.ones(num_nodes))
 
         # Setup partials
@@ -68,6 +71,7 @@ def configure_io(self):
         self.declare_partials(of='t', wrt='t_initial', rows=rs, cols=cs, val=1.0)
         self.declare_partials(of='t', wrt='t_duration', rows=rs, cols=cs, val=dtime_dduration)
         self.declare_partials(of='t_phase', wrt='t_duration', rows=rs, cols=cs, val=dtime_dduration)
+        self.declare_partials(of='t_final',  wrt=['t_initial', 't_duration'], val=1.0)
         self.declare_partials(of='dt_dstau', wrt='t_duration', rows=rs, cols=cs, val=0.5 * node_dptau_dstau)
 
     def compute(self, inputs, outputs):
@@ -90,3 +94,4 @@ def compute(self, inputs, outputs):
         outputs['t'][:] = t_initial + 0.5 * (node_ptau + 1) * t_duration
         outputs['t_phase'][:] = 0.5 * (node_ptau + 1) * t_duration
         outputs['dt_dstau'][:] = 0.5 * t_duration * node_dptau_dstau
+        outputs['t_final'] = t_initial + t_duration
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index d60f2912b..adb64ace2 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -220,7 +220,8 @@ def configure_time(self, phase):
                               flat_src_indices=True if dynamic else None)
 
         for name, tgts in [('t_initial', time_options['t_initial_targets']),
-                           ('t_duration', time_options['t_duration_targets'])]:
+                           ('t_duration', time_options['t_duration_targets']),
+                           ('t_final', time_options['t_final_targets'])]:
 
             targets = _get_targets_metadata(ode, name, user_targets=tgts)
             for t, meta in targets.items():
diff --git a/dymos/transcriptions/grid_data.py b/dymos/transcriptions/grid_data.py
index eb668c4e9..c6da1dd64 100644
--- a/dymos/transcriptions/grid_data.py
+++ b/dymos/transcriptions/grid_data.py
@@ -491,6 +491,16 @@ def __eq__(self, other):
         else:
             return False
 
+    def __str__(self):
+        s = (f'{self.__class__.__name__}(num_seg={self.num_segments}, '
+             f'num_nodes={self.num_nodes}, '
+             f'compressed={self.compressed})\n')
+        s += f'  order = {self.transcription_order}\n'
+        s += f'  segment_ends = {self.segment_ends}\n'
+        s += '  num_nodes_per_segment:\n'
+        s += '\n'.join([f'    {subset_name} = {nnps}' for subset_name, nnps in sorted(self.subset_num_nodes_per_segment.items())])
+        return s
+
     def is_aligned_with(self, other, tol=1.0E-12):
         """
         Check that the segment distribution in GridData object `other` matches that of this GridData object.
diff --git a/dymos/transcriptions/pseudospectral/birkhoff.py b/dymos/transcriptions/pseudospectral/birkhoff.py
index 06a83c9bd..ea8ecaa52 100644
--- a/dymos/transcriptions/pseudospectral/birkhoff.py
+++ b/dymos/transcriptions/pseudospectral/birkhoff.py
@@ -108,11 +108,10 @@ def configure_time(self, phase):
                 src_idxs = self.grid_data.subset_node_indices['all']
                 phase.connect(name, [f'ode_all.{t}' for t in targets], src_indices=src_idxs,
                               flat_src_indices=True)
-                src_idxs = om.slicer[[0, -1], ...]
-                phase.connect(name, [f'boundary_vals.{t}' for t in targets], src_indices=src_idxs)
 
         for name, targets in [('t_initial', options['t_initial_targets']),
-                              ('t_duration', options['t_duration_targets'])]:
+                              ('t_duration', options['t_duration_targets']),
+                              ('t_final', options['t_final_targets'])]:
             for t in targets:
                 shape = ode_inputs[t]['shape']
 
diff --git a/dymos/transcriptions/pseudospectral/components/__init__.py b/dymos/transcriptions/pseudospectral/components/__init__.py
index 9b9e01c09..a6f8bbc5d 100644
--- a/dymos/transcriptions/pseudospectral/components/__init__.py
+++ b/dymos/transcriptions/pseudospectral/components/__init__.py
@@ -3,6 +3,8 @@
 from .state_independents import StateIndependentsComp
 from .control_endpoint_defect_comp import ControlEndpointDefectComp
 from .gauss_lobatto_interleave_comp import GaussLobattoInterleaveComp
-from .birkhoff_collocation_comp import BirkhoffCollocationComp
+from .birkhoff_defect_comp import BirkhoffDefectComp
 from .birkhoff_iter_group import BirkhoffIterGroup
 from .birkhoff_boundary_group import BirkhoffBoundaryGroup
+from .radau_iter_group import RadauIterGroup
+from .radau_boundary_group import RadauBoundaryGroup
diff --git a/dymos/transcriptions/pseudospectral/components/birkhoff_boundary_group.py b/dymos/transcriptions/pseudospectral/components/birkhoff_boundary_group.py
index 26a79fbba..9d3138fbb 100644
--- a/dymos/transcriptions/pseudospectral/components/birkhoff_boundary_group.py
+++ b/dymos/transcriptions/pseudospectral/components/birkhoff_boundary_group.py
@@ -80,7 +80,7 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
 
 class BirkhoffBoundaryGroup(om.Group):
     """
-    Class definition for the BirkhoffBoundaryEvalGroup.
+    Class definition for the BirkhoffBoundaryGroup.
 
     This group accepts values for initial and final times, states, controls, and parameters
     and evaluates the ODE with those in order to compute the boundary values and
diff --git a/dymos/transcriptions/pseudospectral/components/birkhoff_collocation_comp.py b/dymos/transcriptions/pseudospectral/components/birkhoff_defect_comp.py
similarity index 98%
rename from dymos/transcriptions/pseudospectral/components/birkhoff_collocation_comp.py
rename to dymos/transcriptions/pseudospectral/components/birkhoff_defect_comp.py
index f40ca6aa5..32066e678 100644
--- a/dymos/transcriptions/pseudospectral/components/birkhoff_collocation_comp.py
+++ b/dymos/transcriptions/pseudospectral/components/birkhoff_defect_comp.py
@@ -12,11 +12,11 @@
 from dymos.utils.birkhoff import birkhoff_matrix
 
 
-class BirkhoffCollocationComp(om.ExplicitComponent):
+class BirkhoffDefectComp(om.ExplicitComponent):
     """
-    Class definition for the BirkhoffCollocationComp.
+    Class definition for the BirkhoffDefectComp.
 
-    BirkhoffCollocationComp computes the generalized defects of a segment for implicit collocation.
+    BirkhoffDefectComp computes the generalized defects of a segment for implicit collocation.
     There are four defects to be evaluated; state, state rate, initial state, and final state.
     The state defect is the difference between the state value design variables and the state rate
     design variables.
diff --git a/dymos/transcriptions/pseudospectral/components/birkhoff_iter_group.py b/dymos/transcriptions/pseudospectral/components/birkhoff_iter_group.py
index d26389c79..79179de5b 100644
--- a/dymos/transcriptions/pseudospectral/components/birkhoff_iter_group.py
+++ b/dymos/transcriptions/pseudospectral/components/birkhoff_iter_group.py
@@ -1,8 +1,8 @@
 import numpy as np
 import openmdao.api as om
 
-from .birkhoff_collocation_comp import BirkhoffCollocationComp
-from .birkhoff_state_resid_comp import BirkhoffStateResidComp
+from .birkhoff_defect_comp import BirkhoffDefectComp
+from .input_resids_comp import InputResidsComp
 
 from ...grid_data import GridData
 from ....phase.options import TimeOptionsDictionary
@@ -53,13 +53,13 @@ def setup(self):
         self.add_subsystem('ode_all', subsys=ode_class(num_nodes=nn, **ode_init_kwargs))
 
         self.add_subsystem('collocation_comp',
-                           subsys=BirkhoffCollocationComp(grid_data=gd,
+                           subsys=BirkhoffDefectComp(grid_data=gd,
                                                           state_options=state_options,
                                                           time_units=time_options['units']),
                            promotes_inputs=['*'], promotes_outputs=['*'])
 
         if any([opts['solve_segments'] in ('forward', 'backward') for opts in state_options.values()]):
-            self.add_subsystem('states_balance_comp', subsys=BirkhoffStateResidComp(),
+            self.add_subsystem('states_balance_comp', subsys=InputResidsComp(),
                                promotes_inputs=['*'], promotes_outputs=['*'])
 
     def _configure_desvars(self, name, options):
@@ -208,20 +208,20 @@ def configure_io(self, phase):
                                                shape=(nn,) + shape,
                                                units=units)
 
-                states_balance_comp.add_residual_from_input(f'state_defects:{name}',
-                                                            shape=(nn+ns,) + shape,
-                                                            units=units)
+                states_balance_comp.add_input(f'state_defects:{name}',
+                                              shape=(nn+ns,) + shape,
+                                              units=units)
 
                 if ns > 1:
-                    states_balance_comp.add_residual_from_input(f'state_cnty_defects:{name}',
-                                                                shape=(ns - 1,) + shape,
-                                                                units=units)
+                    states_balance_comp.add_input(f'state_cnty_defects:{name}',
+                                                  shape=(ns - 1,) + shape,
+                                                  units=units)
 
             if f'state_rates:{name}' in self._implicit_outputs:
                 states_balance_comp.add_output(f'state_rates:{name}', shape=(nn,) + shape, units=units)
-                states_balance_comp.add_residual_from_input(f'state_rate_defects:{name}',
-                                                            shape=(nn,) + shape,
-                                                            units=units)
+                states_balance_comp.add_input(f'state_rate_defects:{name}',
+                                              shape=(nn,) + shape,
+                                              units=units)
 
             if f'initial_states:{name}' in self._implicit_outputs:
                 states_balance_comp.add_output(f'initial_states:{name}', shape=(1,) + shape, units=units)
diff --git a/dymos/transcriptions/pseudospectral/components/birkhoff_state_resid_comp.py b/dymos/transcriptions/pseudospectral/components/birkhoff_state_resid_comp.py
deleted file mode 100644
index 09cde3b36..000000000
--- a/dymos/transcriptions/pseudospectral/components/birkhoff_state_resid_comp.py
+++ /dev/null
@@ -1,72 +0,0 @@
-import numpy as np
-
-import openmdao.api as om
-
-from openmdao.utils.general_utils import ensure_compatible
-
-from dymos._options import options as dymos_options
-from dymos.utils.misc import om_version
-
-
-class BirkhoffStateResidComp(om.ImplicitComponent):
-    """
-    Class definition for the BirkhoffStateResidComp.
-
-    Generates the residuals for any states that are solved for implicitly.
-
-    Parameters
-    ----------
-    **kwargs : dict
-        Dictionary of optional arguments.
-    """
-
-    def __init__(self, **kwargs):
-        super().__init__(**kwargs)
-        self._no_check_partials = not dymos_options['include_check_partials']
-        self._io_pairs = []
-
-    def add_residual_from_input(self, name, **kwargs):
-        """
-        Adds a residual whose value is given by resid_input.
-
-        Parameters
-        ----------
-        name : str
-            The name of the input providing the residuals for the given output.
-        **kwargs : dict
-            Additional keyword arguments for add_input and add_residual.
-        """
-        val = kwargs['val'] if 'val' in kwargs else 1.0
-        shape = kwargs['shape'] if 'shape' in kwargs else None
-        val, shape = ensure_compatible(name, value=val, shape=shape)
-        resid_name = 'resid_' + name
-
-        self._io_pairs.append((resid_name, name))
-
-        size = np.prod(shape, dtype=int)
-        ar = np.arange(size, dtype=int)
-
-        self.add_input(name, **kwargs)
-        self.add_residual(resid_name, **kwargs)
-
-        if om_version()[0] > (3, 31, 1):
-            self.declare_partials(of=resid_name, wrt=name, rows=ar, cols=ar, val=1.0)
-        else:
-            self.declare_partials(of='*', wrt='*', method='fd')
-
-    def apply_nonlinear(self, inputs, outputs, residuals):
-        """
-        Compute residuals given inputs and outputs.
-
-        The model is assumed to be in an unscaled state.
-
-        Parameters
-        ----------
-        inputs : Vector
-            Unscaled, dimensional input variables read via inputs[key].
-        outputs : Vector
-            Unscaled, dimensional output variables read via outputs[key].
-        residuals : Vector
-            Unscaled, dimensional residuals written to via residuals[key].
-        """
-        residuals.set_val(inputs.asarray())
diff --git a/dymos/transcriptions/pseudospectral/components/input_resids_comp.py b/dymos/transcriptions/pseudospectral/components/input_resids_comp.py
new file mode 100644
index 000000000..035be2564
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/input_resids_comp.py
@@ -0,0 +1,117 @@
+"""InputResidsComp provides a simple implicit component with minimal boilerplate."""
+
+import numpy as np
+
+from openmdao.core.implicitcomponent import ImplicitComponent
+from openmdao.utils.general_utils import ensure_compatible
+
+
+# TODO: Remove this component when the minimum supported OpenMDAO version is >3.35.0
+
+class InputResidsComp(ImplicitComponent):
+    """
+    Class definition for the InputResidsComp.
+
+    Uses all inputs as residuals while allowing individual outputs that are not necessarily
+    associated with a specific residual.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+
+    Attributes
+    ----------
+    _refs : dict
+        Residual ref values that are cached during calls to the overloaded add_input method.
+    """
+
+    def __init__(self, **kwargs):
+        """
+        Initialize the InputResidsComp.
+
+        Parameters
+        ----------
+        **kwargs : dict
+            Keyword arguments passed to the __init__ method of ImplicitComponent
+        """
+        self._refs = {}
+        super().__init__(**kwargs)
+
+    def add_input(self, name, val=1.0, shape=None, units=None, desc='', tags=None,
+                  shape_by_conn=False, copy_shape=None, compute_shape=None, distributed=None,
+                  ref=None):
+        """
+        Add an input to be used as a residual.
+
+        Parameters
+        ----------
+        name : str
+            Name of the variable in this component's namespace.
+        val : float or list or tuple or ndarray or Iterable
+            The initial value of the variable being added in user-defined units.
+            Default is 1.0.
+        shape : int or tuple or list or None
+            Shape of this variable, only required if val is not an array. Default is None.
+        units : str or None
+            Units in which this input variable will be provided to the component
+            during execution. Default is None, which means it is unitless.
+        desc : str
+            Description of the variable.
+        tags : str or list of strs
+            User defined tags that can be used to filter what gets listed when calling
+            list_inputs and list_outputs.
+        shape_by_conn : bool
+            If True, shape this input to match its connected output.
+        copy_shape : str or None
+            If a str, that str is the name of a variable. Shape this input to match that of
+            the named variable.
+        compute_shape : function
+            A function taking a dict arg containing names and shapes of this component's outputs
+            and returning the shape of this input.
+        distributed : bool
+            If True, this variable is a distributed variable, so it can have different sizes/values
+            across MPI processes.
+        ref : float or ndarray or None
+            Scaling parameter. The value in the user-defined units of this residual
+            when the scaled value is 1. Default is 1.
+        """
+        self._refs[name] = ref
+        super().add_input(name, val=val, shape=shape, units=units, desc=desc, tags=tags,
+                          shape_by_conn=shape_by_conn, copy_shape=copy_shape,
+                          compute_shape=compute_shape, distributed=distributed)
+
+        val, shape = ensure_compatible(name, val, shape)
+
+        self.add_residual(f'resid_{name}', shape=shape, units=units,
+                          ref=self._refs[name])
+
+    def setup_partials(self):
+        """
+        Delay calls to declare_partials for the component.
+
+        This method is used because input/residual sizes
+        may not be known until final setup.
+        """
+        for name in self._var_rel_names['input']:
+            resid_name = 'resid_' + name
+            size = self._var_rel2meta[name]['size']
+            ar = np.arange(size, dtype=int)
+            self.declare_partials(of=resid_name, wrt=name, rows=ar, cols=ar, val=1.0)
+
+    def apply_nonlinear(self, inputs, outputs, residuals):
+        """
+        Compute residuals given inputs and outputs.
+
+        The model is assumed to be in an unscaled state.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        outputs : Vector
+            Unscaled, dimensional output variables read via outputs[key].
+        residuals : Vector
+            Unscaled, dimensional residuals written to via residuals[key].
+        """
+        residuals.set_val(inputs.asarray())
diff --git a/dymos/transcriptions/pseudospectral/components/radau_boundary_group.py b/dymos/transcriptions/pseudospectral/components/radau_boundary_group.py
new file mode 100644
index 000000000..e84c14606
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/radau_boundary_group.py
@@ -0,0 +1,149 @@
+import numpy as np
+import openmdao.api as om
+
+from ...grid_data import GridData
+from dymos._options import options as dymos_options
+
+
+class RadauBoundaryMuxComp(om.ExplicitComponent):
+    """
+    Class definition of the RadauBoundaryMuxComp.
+
+    This component takes the initial and final values of states
+    and muxes them into a single output.
+
+    For a state of a given `shape` in a phase with `num_seg` segments,
+    the shape of `initial_states:{state_name}` and `final_states:{state_name}`
+    are both `(num_seg,) + shape` and the shape of the resulting
+    `states:{state_name}` is `(2,) + shape`.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional phase arguments.
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._io_names = {}
+        self._no_check_partials = not dymos_options['include_check_partials']
+
+    def configure_io(self, state_options):
+        """
+        I/O creation is delayed until configure so that we can determine shape and units for the states.
+
+        Parameters
+        ----------
+        state_options : StateOptionsDictionary
+            The phase object to which this transcription instance applies.
+        """
+        self._io_names = {}
+        for state_name, options in state_options.items():
+            shape = options['shape']
+            size = np.prod(shape, dtype=int)
+            units = options['units']
+            self._io_names[state_name] = {'initial': f'initial_states:{state_name}',
+                                          'final': f'final_states:{state_name}',
+                                          'boundary': state_name}
+
+            iname = self._io_names[state_name]['initial']
+            fname = self._io_names[state_name]['final']
+            bname = self._io_names[state_name]['boundary']
+
+            self.add_input(iname, shape=shape, units=units)
+            self.add_input(fname, shape=shape, units=units)
+            self.add_output(bname, shape=(2,) + shape, units=units)
+
+            ar = np.arange(size, dtype=int)
+
+            self.declare_partials(of=bname, wrt=iname, rows=ar, cols=ar, val=1.0)
+            self.declare_partials(of=bname, wrt=fname, rows=ar + size, cols=ar, val=1.0)
+
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        """
+        Compute component outputs.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        outputs : Vector
+            Unscaled, dimensional output variables read via outputs[key].
+        discrete_inputs : dict or None
+            If not None, dict containing discrete input values.
+        discrete_outputs : dict or None
+            If not None, dict containing discrete output values.
+        """
+        for state_name, io_names in self._io_names.items():
+            outputs[io_names['boundary']][0] = inputs[io_names['initial']]
+            outputs[io_names['boundary']][1] = inputs[io_names['final']]
+
+
+class RadauBoundaryGroup(om.Group):
+    """
+    Class definition for the RadauBoundaryGroup.
+
+    This group accepts values for initial and final times, states, controls, and parameters
+    and evaluates the ODE with those in order to compute the boundary values and
+    objectives.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+    """
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._no_check_partials = not dymos_options['include_check_partials']
+
+    def initialize(self):
+        """
+        Declare group options.
+        """
+        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={}, 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.')
+
+    def setup(self):
+        """
+        Define the structure of the control group.
+        """
+        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=RadauBoundaryMuxComp(),
+                           promotes_inputs=['*'], promotes_outputs=['*'])
+
+        self.add_subsystem('boundary_ode', subsys=ode_class(num_nodes=2, **ode_init_kwargs),
+                           promotes_inputs=['*'], promotes_outputs=['*'])
+
+        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=['*'])
+
+    def configure_io(self, phase):
+        """
+        I/O creation is delayed until configure so that we can determine shape and units for the states.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        self._get_subsystem('boundary_mux').configure_io(state_options=phase.state_options)
+
+        for state_name, options in phase.state_options.items():
+            for tgt in options['targets']:
+                self.promotes('boundary_ode', inputs=[(tgt, state_name)])
diff --git a/dymos/transcriptions/pseudospectral/components/radau_defect_comp.py b/dymos/transcriptions/pseudospectral/components/radau_defect_comp.py
new file mode 100644
index 000000000..99e93bd7f
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/radau_defect_comp.py
@@ -0,0 +1,318 @@
+import numpy as np
+import openmdao.api as om
+
+import scipy.sparse as sp
+
+from dymos._options import options as dymos_options
+from dymos.transcriptions.grid_data import GridData
+from dymos.utils.misc import get_rate_units
+
+
+class RadauDefectComp(om.ExplicitComponent):
+    """
+    Class definiton for the RadauDefectComp.
+
+    RadauDefectComp computes the generalized defect of a segment for implicit collocation.
+    The defect is the interpolated state derivative at the collocation nodes minus
+    the computed state derivative at the collocation nodes.
+
+    In addition to specifying the initial and final state values at the endpoints via the
+    polynomial-interpolated `states:{state_name}`, the initial and final state values are
+    also specified via `initial_states:{state_name}` and `final_states:{state_name}`.
+    These difference in these initial and final values are then computed as
+    `initial_state_defects:{state_name}` and `final_state_defects:{state_name}`.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._no_check_partials = not dymos_options['include_check_partials']
+
+    def initialize(self):
+        """Declare component options."""
+        self.options.declare(
+            'grid_data', types=GridData,
+            desc='Container object for grid info')
+
+        self.options.declare(
+            'state_options', types=dict,
+            desc='Dictionary of state names/options for the phase')
+
+        self.options.declare(
+            'time_units', default=None, allow_none=True, types=str,
+            desc='Units of time')
+
+    def configure_io(self, phase):
+        """
+        I/O creation is delayed until configure so we can determine shape and units.
+
+        Parameters
+        ----------
+        phase : Phase
+            The phase in which this component exists.
+        """
+        gd: GridData = self.options['grid_data']
+        num_segs: int = gd.num_segments
+        num_nodes: int = gd.subset_num_nodes['all']
+        num_col_nodes: int = gd.subset_num_nodes['col']
+        time_units: str = self.options['time_units']
+        state_options = self.options['state_options']
+
+        # The radau differentiation matrix
+        _, self._D = self.options['grid_data'].phase_lagrange_matrices('state_disc',
+                                                                       'col',
+                                                                       sparse=True)
+
+        self.add_input('dt_dstau', units=time_units, shape=(num_col_nodes,))
+
+        self.var_names = var_names = {}
+        for state_name in state_options:
+            var_names[state_name] = {
+                'initial_val': f'initial_states:{state_name}',
+                'final_val': f'final_states:{state_name}',
+                'val': f'states:{state_name}',
+                'f_ode': f'f_ode:{state_name}',
+                'rate_defect': f'state_rate_defects:{state_name}',
+                'cnty_defect': f'state_cnty_defects:{state_name}',
+                'initial_defect': f'initial_state_defects:{state_name}',
+                'final_defect': f'final_state_defects:{state_name}'
+            }
+
+        for state_name, options in state_options.items():
+            shape = options['shape']
+            units = options['units']
+
+            rate_units = get_rate_units(units, time_units)
+
+            var_names = self.var_names[state_name]
+
+            self.add_input(name=var_names['initial_val'],
+                           shape=(1,) + shape,
+                           units=units,
+                           desc='Initial value of the state at the start of the phase.')
+
+            self.add_input(name=var_names['final_val'],
+                           shape=(1,) + shape,
+                           units=units,
+                           desc='Final value of the state at the end of the phase.')
+
+            self.add_input(var_names['val'],
+                           shape=(num_nodes,) + shape,
+                           units=units,
+                           desc='state value at all nodes within the phase')
+
+            self.add_input(
+                name=var_names['f_ode'],
+                shape=(num_col_nodes,) + shape,
+                desc=f'Computed derivative of state {state_name} at the collocation nodes',
+                units=rate_units)
+
+            self.add_output(
+                name=var_names['initial_defect'],
+                shape=(1,) + shape,
+                desc=f'Initial value defect of state {state_name}',
+                units=units)
+
+            self.add_output(
+                name=var_names['final_defect'],
+                shape=(1,) + shape,
+                desc=f'Final value defect of state {state_name}',
+                units=units)
+
+            self.add_output(
+                name=var_names['rate_defect'],
+                shape=(num_col_nodes,) + shape,
+                desc=f'Interior defects of state {state_name}',
+                units=units)
+
+            if gd.num_segments > 1 and not gd.compressed:
+                self.add_output(
+                    name=var_names['cnty_defect'],
+                    shape=(num_segs - 1,) + shape,
+                    desc=f'Segment boundary defect of state {state_name}',
+                    units=units)
+
+            if 'defect_ref' in options and options['defect_ref'] is not None:
+                defect_ref = options['defect_ref']
+            elif 'defect_scaler' in options and options['defect_scaler'] is not None:
+                defect_ref = np.divide(1.0, options['defect_scaler'])
+            elif 'ref' in options and options['ref'] is not None:
+                defect_ref = options['ref']
+            elif 'scaler' in options and options['scaler'] is not None:
+                defect_ref = np.divide(1.0, options['scaler'])
+            else:
+                defect_ref = 1.0
+
+            if np.isscalar(defect_ref):
+                defect_ref = defect_ref * np.ones(shape)
+
+            if defect_ref.shape != shape:
+                raise ValueError('array-valued scaler/ref/defect_ref must be the same shape as the state')
+
+            rate_defect_ref = np.tile(defect_ref, (num_col_nodes, 1))
+
+            if not options['solve_segments']:
+                self.add_constraint(name=var_names['rate_defect'],
+                                    equals=0.0,
+                                    ref=rate_defect_ref)
+
+                self.add_constraint(name=var_names['initial_defect'],
+                                    equals=0.0,
+                                    ref=defect_ref)
+
+                self.add_constraint(name=var_names['final_defect'],
+                                    equals=0.0,
+                                    ref=defect_ref)
+
+                if gd.num_segments > 1 and not gd.compressed:
+                    cnty_defect_ref = np.tile(np.asarray(defect_ref), (num_segs - 1, 1))
+                    self.add_constraint(name=var_names['cnty_defect'],
+                                        equals=0.0,
+                                        ref=cnty_defect_ref)
+
+    def setup_partials(self):
+        gd: GridData = self.options['grid_data']
+        num_segs: int = gd.num_segments
+        num_nodes: int = gd.subset_num_nodes['all']
+        num_col_nodes: int = gd.subset_num_nodes['col']
+        state_options = self.options['state_options']
+
+        num_col_nodes = self.options['grid_data'].subset_num_nodes['col']
+        state_options = self.options['state_options']
+
+        for state_name, options in state_options.items():
+            shape = options['shape']
+            size = np.prod(shape)
+
+            r = np.arange(num_col_nodes * size)
+
+            var_names = self.var_names[state_name]
+
+            self.declare_partials(of=var_names['rate_defect'],
+                                  wrt=var_names['f_ode'],
+                                  rows=r, cols=r, val=-1.0)
+
+            c = np.repeat(np.arange(num_col_nodes), size)
+            self.declare_partials(of=var_names['rate_defect'],
+                                  wrt='dt_dstau',
+                                  rows=r, cols=c)
+
+            # The state rate defects wrt the state values at the discretization nodes
+            # are given by the differentiation matrix.
+            sparse_D_of_size = sp.kron(sp.csr_matrix(self._D), sp.eye(size), format='coo')
+            r, c = sparse_D_of_size.nonzero()
+
+            self.declare_partials(of=var_names['rate_defect'],
+                                  wrt=var_names['val'],
+                                  rows=r,
+                                  cols=c,
+                                  val=sparse_D_of_size.data)
+
+            # The initial value defect is just an identity matrix at the "top left" corner of the jacobian.
+            ar_size = np.arange(size, dtype=int)
+            self.declare_partials(of=var_names['initial_defect'],
+                                  wrt=var_names['val'],
+                                  rows=ar_size, cols=ar_size, val=-1.0)
+
+            self.declare_partials(of=var_names['initial_defect'],
+                                  wrt=var_names['initial_val'],
+                                  rows=ar_size, cols=ar_size, val=1.0)
+
+            # The final value defect is an identity matrix at the "bottom right" corner of the jacobian.
+            # r = np.arange(size, dtype=int)
+            # c = np.arange(num_nodes - size, num_nodes, dtype=int)
+            row_vec_end_1 = np.zeros((1, num_nodes))
+            row_vec_end_1[:, -1] = -1.0
+            pattern = sp.kron(row_vec_end_1, sp.eye(size), format='coo')
+            r, c = pattern.nonzero()
+            self.declare_partials(of=var_names['final_defect'],
+                                  wrt=var_names['val'],
+                                  rows=r, cols=c, val=-1.0)
+
+            self.declare_partials(of=var_names['final_defect'],
+                                  wrt=var_names['final_val'],
+                                  rows=ar_size, cols=ar_size, val=1.0)
+
+            if gd.num_segments > 1 and not gd.compressed:
+                idxs_se: int = gd.subset_node_indices['segment_ends']
+                seg_end_pattern = np.zeros((num_segs - 1, num_nodes), dtype=int)
+                for i_row in range(num_segs - 1):
+                    end_idx = idxs_se[1:-1][2 * i_row]
+                    start_idx = idxs_se[1:-1][2 * i_row + 1]
+                    seg_end_pattern[i_row, end_idx] = -1
+                    seg_end_pattern[i_row, start_idx] = 1
+
+                pattern = np.kron(seg_end_pattern, np.eye(size, dtype=int))
+
+                r, c = pattern.nonzero()
+
+                self.declare_partials(of=var_names['cnty_defect'],
+                                      wrt=var_names['val'], rows=r, cols=c, val=pattern[r, c])
+
+    def compute(self, inputs, outputs):
+        """
+        Compute collocation defects.
+
+        Parameters
+        ----------
+        inputs : `Vector`
+            `Vector` containing inputs.
+        outputs : `Vector`
+            `Vector` containing outputs.
+        """
+        gd: GridData = self.options['grid_data']
+        num_disc_nodes: int = gd.subset_num_nodes['state_disc']
+        num_col_nodes: int = gd.subset_num_nodes['col']
+        idxs_se: int = gd.subset_node_indices['segment_ends']
+
+        state_options = self.options['state_options']
+        dt_dstau: np.ndarray = inputs['dt_dstau']
+        D = self._D
+
+        for state_name, state_options in state_options.items():
+            shape = state_options['shape']
+            size = np.prod(shape)
+            var_names = self.var_names[state_name]
+
+            f_ode = inputs[var_names['f_ode']]
+            x = inputs[var_names['val']]
+            x_0 = inputs[var_names['initial_val']]
+            x_f = inputs[var_names['final_val']]
+
+            # The defect is computed as
+            # defect = D @ x - f_ode * dt_dstau  # noqa: ERA001
+            # But scipy.sparse only handles 2D matrices, so we need to force x to be 2D
+            # and then change the product back to the proper shape.
+
+            x_flat = np.reshape(x, newshape=(num_disc_nodes, size))
+            f_approx = np.reshape(D.dot(x_flat), newshape=(num_col_nodes,) + shape)
+
+            outputs[var_names['rate_defect']] = f_approx - (f_ode.T * dt_dstau).T
+            outputs[var_names['initial_defect']] = x_0 - x[0, ...]
+            outputs[var_names['final_defect']] = x_f - x[-1, ...]
+
+            if gd.num_segments > 1 and not gd.compressed:
+                outputs[var_names['cnty_defect']] = x[idxs_se[2::2], ...] - x[idxs_se[1:-2:2], ...]
+
+    def compute_partials(self, inputs, partials):
+        """
+        Compute sub-jacobian parts. The model is assumed to be in an unscaled state.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        partials : Jacobian
+            Subjac components written to partials[output_name, input_name].
+        """
+        dt_dstau = inputs['dt_dstau']
+        for state_name, options in self.options['state_options'].items():
+            size = np.prod(options['shape'])
+            var_names = self.var_names[state_name]
+            f_ode = inputs[var_names['f_ode']]
+
+            partials[var_names['rate_defect'], var_names['f_ode']] = -np.repeat(dt_dstau, size)
+            partials[var_names['rate_defect'], 'dt_dstau'] = -f_ode.ravel()
diff --git a/dymos/transcriptions/pseudospectral/components/radau_iter_group.py b/dymos/transcriptions/pseudospectral/components/radau_iter_group.py
new file mode 100644
index 000000000..a969cbd78
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/radau_iter_group.py
@@ -0,0 +1,356 @@
+import numpy as np
+import openmdao.api as om
+
+from .input_resids_comp import InputResidsComp
+from .radau_defect_comp import RadauDefectComp
+
+from dymos.transcriptions.grid_data import GridData
+from dymos.phase.options import TimeOptionsDictionary
+
+
+class RadauIterGroup(om.Group):
+    """
+    Class definition for the RadauIterGroup.
+
+    This group allows for iteration of the state variables and initial _or_ final value of the state
+    depending on the direction of the solve.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._implicit_outputs = set()
+
+    def initialize(self):
+        """Declare group options."""
+        self.options.declare('state_options', types=dict,
+                             desc='Dictionary of options for the states.')
+        self.options.declare('time_options', types=TimeOptionsDictionary,
+                             desc='Options for time in 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')
+
+    def setup(self):
+        """
+        Define the structure of the RadauIterGroup.
+        """
+        gd = self.options['grid_data']
+        nn = gd.subset_num_nodes['all']
+        state_options = self.options['state_options']
+        time_options = self.options['time_options']
+        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))
+
+        self.add_subsystem('defects',
+                           subsys=RadauDefectComp(grid_data=gd,
+                                                  state_options=state_options,
+                                                  time_units=time_options['units']),
+                           promotes_inputs=['*'], promotes_outputs=['*'])
+
+        if any([opts['solve_segments'] in ('forward', 'backward') for opts in state_options.values()]):
+            self.add_subsystem('states_resids_comp', subsys=InputResidsComp(),
+                               promotes_inputs=['*'], promotes_outputs=['*'])
+
+    def _configure_desvars(self, name, options):
+        state_name = f'states:{name}'
+        initial_state_name = f'initial_states:{name}'
+        final_state_name = f'final_states:{name}'
+        state_rate_name = f'state_rates:{name}'
+
+        solve_segs = options['solve_segments']
+        opt = options['opt']
+
+        num_nodes = self.options['grid_data'].subset_num_nodes['col']
+
+        ib = (None, None) if options['initial_bounds'] is None else options['initial_bounds']
+        fb = (None, None) if options['final_bounds'] is None else options['final_bounds']
+        lower = options['lower']
+        upper = options['upper']
+        scaler = options['scaler']
+        adder = options['adder']
+        ref0 = options['ref0']
+        ref = options['ref']
+        fix_initial = options['fix_initial']
+        fix_final = options['fix_final']
+        input_initial = options['input_initial']
+        input_final = options['input_final']
+        shape = options['shape']
+
+        if solve_segs == 'forward' and fix_final:
+            raise ValueError(f"Option fix_final on state {name} may not "
+                             f"be used with `solve_segments='forward'`.\n Use "
+                             f"a boundary constraint to constrain the final "
+                             f"state value instead.")
+        elif solve_segs == 'backward' and fix_initial:
+            raise ValueError(f"Option fix_final on state {name} may not "
+                             f"be used with `solve_segments='forward'`.\n Use "
+                             f"a boundary constraint to constrain the initial "
+                             f"state value instead.")
+
+        if not np.isscalar(ref0) and ref0 is not None:
+            ref0 = np.asarray(ref0)
+            if ref0.shape == shape:
+                ref0_state = np.tile(ref0.flatten(), num_nodes)
+            else:
+                raise ValueError('array-valued scaler/ref must length equal to state-size')
+        else:
+            ref0_state = ref0
+        if not np.isscalar(ref) and ref is not None:
+            ref = np.asarray(ref)
+            if ref.shape == shape:
+                ref_state = np.tile(ref.flatten(), num_nodes)
+            else:
+                raise ValueError('array-valued scaler/ref must length equal to state-size')
+        else:
+            ref_state = ref
+
+        free_vars = {state_name, initial_state_name, final_state_name}
+
+        if solve_segs == 'forward':
+            implicit_outputs = {state_name, final_state_name}
+        elif solve_segs == 'backward':
+            implicit_outputs = {state_name, initial_state_name}
+        else:
+            implicit_outputs = set()
+
+        free_vars = free_vars - implicit_outputs
+
+        if fix_initial or input_initial:
+            free_vars = free_vars - {initial_state_name}
+        if fix_final or input_final:
+            free_vars = free_vars - {final_state_name}
+
+        if opt:
+            # Add design variables for the remaining free variables
+            if state_name in free_vars:
+                self.add_design_var(name=state_name,
+                                    lower=lower,
+                                    upper=upper,
+                                    scaler=scaler,
+                                    adder=adder,
+                                    ref0=ref0_state,
+                                    ref=ref_state)
+
+            if state_rate_name in free_vars:
+                self.add_design_var(name=state_rate_name,
+                                    scaler=scaler,
+                                    adder=adder,
+                                    ref0=ref0_state,
+                                    ref=ref_state)
+
+            if initial_state_name in free_vars:
+                self.add_design_var(name=initial_state_name,
+                                    lower=ib[0],
+                                    upper=ib[1],
+                                    scaler=scaler,
+                                    adder=adder,
+                                    ref0=ref0,
+                                    ref=ref)
+
+            if final_state_name in free_vars:
+                self.add_design_var(name=final_state_name,
+                                    lower=fb[0],
+                                    upper=fb[1],
+                                    scaler=scaler,
+                                    adder=adder,
+                                    ref0=ref0,
+                                    ref=ref)
+
+        return implicit_outputs
+
+    def configure_io(self, phase):
+        """
+        I/O creation is delayed until configure so that we can determine shape and units for the states.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        defect_comp = self._get_subsystem('defects')
+        defect_comp.configure_io(phase)
+
+        gd = self.options['grid_data']
+        nn = gd.subset_num_nodes['all']
+        nin = gd.subset_num_nodes['state_input']
+        ncn = gd.subset_num_nodes['col']
+        ns = gd.num_segments
+        state_src_idxs_input_to_disc = gd.input_maps['state_input_to_disc']
+        state_src_idxs_input_to_all = state_src_idxs_input_to_disc
+
+        col_idxs = gd.subset_node_indices['col']
+
+        state_options = self.options['state_options']
+        states_resids_comp = self._get_subsystem('states_resids_comp')
+
+        self.promotes('defects', inputs=('dt_dstau',),
+                      src_indices=om.slicer[col_idxs, ...],
+                      src_shape=(nn,))
+
+        for name, options in state_options.items():
+            units = options['units']
+            rate_source = options['rate_source']
+            shape = options['shape']
+
+            for tgt in options['targets']:
+                self.promotes('ode_all', [(tgt, f'states:{name}')],
+                              src_indices=om.slicer[state_src_idxs_input_to_all, ...])
+
+            self.set_input_defaults(f'states:{name}', val=1.0, units=units, src_shape=(nin,) + shape)
+
+            self.promotes('defects', inputs=(f'states:{name}',),
+                          src_indices=om.slicer[state_src_idxs_input_to_all, ...],)
+
+            self._implicit_outputs = self._configure_desvars(name, options)
+
+            if f'states:{name}' in self._implicit_outputs:
+
+                states_resids_comp.add_input(f'initial_state_defects:{name}', shape=(1,) + shape, units=units)
+                states_resids_comp.add_input(f'final_state_defects:{name}', shape=(1,) + shape, units=units)
+                states_resids_comp.add_input(f'state_rate_defects:{name}', shape=(ncn,) + shape, units=units)
+
+                if ns > 1 and not gd.compressed:
+                    states_resids_comp.add_input(f'state_cnty_defects:{name}',
+                                                 shape=(ns - 1,) + shape,
+                                                 units=units)
+                    # For compressed transcription, resids does not provide values at overlapping
+                    # segment boundaries.
+                    states_resids_comp.add_output(f'states:{name}',
+                                                  shape=(nn,) + shape,
+                                                  units=units)
+                else:
+                    # For compressed transcirption, resids comp provides values at input nodes.
+                    nin = gd.subset_num_nodes['state_input']
+                    states_resids_comp.add_output(f'states:{name}',
+                                                  shape=(nin,) + shape,
+                                                  units=units)
+
+            if f'initial_states:{name}' in self._implicit_outputs:
+                # states_resids_comp.add_input(f'initial_state_defects:{name}', shape=(1,) + shape, units=units)
+                states_resids_comp.add_output(f'initial_states:{name}', shape=(1,) + shape, units=units)
+
+            if f'final_states:{name}' in self._implicit_outputs:
+                # states_resids_comp.add_input(f'final_state_defects:{name}', shape=(1,) + shape, units=units)
+                states_resids_comp.add_output(f'final_states:{name}', shape=(1,) + shape, units=units)
+
+            try:
+                rate_source_var = options['rate_source']
+            except RuntimeError:
+                raise ValueError(f"state '{name}' in phase '{phase.name}' was not given a "
+                                 "rate_source")
+
+            # Note the rate source must be shape-compatible with the state
+            var_type = phase.classify_var(rate_source_var)
+
+            if var_type == 'ode':
+                self.connect(f'ode_all.{rate_source}', f'f_ode:{name}',
+                             src_indices=om.slicer[gd.subset_node_indices['col'], ...])
+
+    def _get_rate_source_path(self, state_name, nodes, phase):
+        """
+        Return the rate source location and indices for a given state name.
+
+        Parameters
+        ----------
+        state_name : str
+            Name of the state.
+        nodes : str
+            One of ['col', 'all'].
+        phase : dymos.Phase
+            Phase object containing the rate source.
+
+        Returns
+        -------
+        str
+            Path to the rate source.
+        ndarray
+            Array of source indices.
+
+        """
+        gd = self.grid_data
+        try:
+            var = phase.state_options[state_name]['rate_source']
+        except RuntimeError:
+            raise ValueError(f"state '{state_name}' in phase '{phase.name}' was not given a "
+                             "rate_source")
+
+        # Note the rate source must be shape-compatible with the state
+        var_type = phase.classify_var(var)
+
+        # Determine the path to the variable
+        if var_type == 't':
+            rate_path = 't'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 't_phase':
+            rate_path = 't_phase'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'state':
+            rate_path = f'states:{var}'
+            # Find the state_input indices which occur at segment endpoints, and repeat them twice
+            state_input_idxs = gd.subset_node_indices['state_input']
+            repeat_idxs = np.ones_like(state_input_idxs)
+            if gd.compressed:
+                segment_end_idxs = gd.subset_node_indices['segment_ends'][1:-1]
+                # Repeat nodes that are on segment bounds (but not the first or last nodes in the phase)
+                nodes_to_repeat = list(set(state_input_idxs).intersection(segment_end_idxs))
+                # Now find these nodes in the state input indices
+                idxs_of_ntr_in_state_inputs = np.where(np.isin(state_input_idxs, nodes_to_repeat))[0]
+                # All state input nodes are used once, but nodes_to_repeat are used twice
+                repeat_idxs[idxs_of_ntr_in_state_inputs] = 2
+            # Now we have a way of mapping the state input indices to all nodes
+            map_input_node_idxs_to_all = np.repeat(np.arange(gd.subset_num_nodes['state_input'],
+                                                             dtype=int), repeats=repeat_idxs)
+            # Now select the subset of nodes we want to use.
+            node_idxs = map_input_node_idxs_to_all[gd.subset_node_indices[nodes]]
+        elif var_type == 'indep_control':
+            rate_path = f'control_values:{var}'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'input_control':
+            rate_path = f'control_values:{var}'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'control_rate':
+            control_name = var[:-5]
+            rate_path = f'control_rates:{control_name}_rate'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'control_rate2':
+            control_name = var[:-6]
+            rate_path = f'control_rates:{control_name}_rate2'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'indep_polynomial_control':
+            rate_path = f'control_values:{var}'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'input_polynomial_control':
+            rate_path = f'control_values:{var}'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'polynomial_control_rate':
+            control_name = var[:-5]
+            rate_path = f'control_rates:{control_name}_rate'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'polynomial_control_rate2':
+            control_name = var[:-6]
+            rate_path = f'control_rates:{control_name}_rate2'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'parameter':
+            rate_path = f'parameter_vals:{var}'
+            dynamic = not phase.parameter_options[var]['static_target']
+            if dynamic:
+                node_idxs = np.zeros(gd.subset_num_nodes[nodes], dtype=int)
+            else:
+                node_idxs = np.zeros(1, dtype=int)
+        else:
+            # Failed to find variable, assume it is in the ODE
+            rate_path = f'ode_all.{var}'
+            node_idxs = gd.subset_node_indices[nodes]
+
+        src_idxs = om.slicer[node_idxs, ...]
+
+        return rate_path, src_idxs
diff --git a/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_collocation_defect.py b/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_collocation_defect.py
index 2df44cfe9..c6b83cd60 100644
--- a/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_collocation_defect.py
+++ b/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_collocation_defect.py
@@ -7,7 +7,7 @@
 from dymos.utils.testing_utils import assert_check_partials
 
 import dymos as dm
-from dymos.transcriptions.pseudospectral.components.birkhoff_collocation_comp import BirkhoffCollocationComp
+from dymos.transcriptions.pseudospectral.components.birkhoff_defect_comp import BirkhoffDefectComp
 from dymos.transcriptions.grid_data import BirkhoffGrid
 
 # Modify class so we can run it standalone.
@@ -23,7 +23,7 @@ def classify_var(self, name):
         return None
 
 
-CollocationComp = CompWrapperConfig(BirkhoffCollocationComp, [_PhaseStub()])
+CollocationComp = CompWrapperConfig(BirkhoffDefectComp, [_PhaseStub()])
 
 
 class TestCollocationComp(unittest.TestCase):
diff --git a/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_iter_group.py b/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_iter_group.py
index 8ffe6070d..7ed817f61 100644
--- a/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_iter_group.py
+++ b/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_iter_group.py
@@ -9,56 +9,13 @@
 from openmdao.utils.testing_utils import use_tempdirs
 
 from dymos.utils.misc import GroupWrapperConfig
+from dymos.utils.testing_utils import PhaseStub, SimpleODE
 from dymos.transcriptions.pseudospectral.components import BirkhoffIterGroup
 from dymos.phase.options import StateOptionsDictionary, TimeOptionsDictionary
 from dymos.transcriptions.grid_data import BirkhoffGrid
 
 
-class _PhaseStub():
-    """
-    A stand-in for the Phase during config_io for testing.
-    It just supports the classify_var method and returns "state", the only value needed for unittests.
-    """
-    def __init__(self):
-        self.nonlinear_solver = None
-        self.linear_solver = None
-
-    def classify_var(self, name):
-        return 'ode'
-
-
-BirkhoffIterGroup = GroupWrapperConfig(BirkhoffIterGroup, [_PhaseStub()])
-
-
-class SimpleODE(om.ExplicitComponent):
-    """
-    A simple ODE from https://math.okstate.edu/people/yqwang/teaching/math4513_fall11/Notes/rungekutta.pdf
-    """
-    def initialize(self):
-        self.options.declare('num_nodes', types=(int,))
-
-    def setup(self):
-        nn = self.options['num_nodes']
-        self.add_input('x', shape=(nn,), units='s**2')
-        self.add_input('t', shape=(nn,), units='s')
-        self.add_input('p', shape=(nn,), units='s**2')
-
-        self.add_output('x_dot', shape=(nn,), units='s')
-
-        ar = np.arange(nn, dtype=int)
-        self.declare_partials(of='x_dot', wrt='x', rows=ar, cols=ar, val=1.0)
-        self.declare_partials(of='x_dot', wrt='t', rows=ar, cols=ar)
-        self.declare_partials(of='x_dot', wrt='p', rows=ar, cols=ar, val=1.0)
-
-    def compute(self, inputs, outputs):
-        x = inputs['x']
-        t = inputs['t']
-        p = inputs['p']
-        outputs['x_dot'] = x - t**2 + p
-
-    def compute_partials(self, inputs, partials):
-        t = inputs['t']
-        partials['x_dot', 't'] = -2*t
+BirkhoffIterGroup = GroupWrapperConfig(BirkhoffIterGroup, [PhaseStub()])
 
 
 @use_tempdirs
diff --git a/dymos/transcriptions/pseudospectral/components/test/test_collocation_defect_solver.py b/dymos/transcriptions/pseudospectral/components/test/test_collocation_defect_solver.py
index f0c41f3c2..9d715c494 100644
--- a/dymos/transcriptions/pseudospectral/components/test/test_collocation_defect_solver.py
+++ b/dymos/transcriptions/pseudospectral/components/test/test_collocation_defect_solver.py
@@ -13,246 +13,6 @@
 from dymos.examples.brachistochrone.test.ex_brachistochrone import brachistochrone_min_time as brach
 
 
-@use_tempdirs
-class TestCollocationBalanceIndex(unittest.TestCase):
-    """
-    Test that the indices used in the StateIndependentsComp are as expected.
-    """
-
-    def setUp(self):
-        dm.options['include_check_partials'] = True
-
-    def tearDown(self):
-        dm.options['include_check_partials'] = False
-
-    def test_3_lgl(self):
-        """
-        Test one 3rd order LGL segment indices
-
-        All nodes are:   [0, 1, 2]
-        Input nodes:     [^     ^]
-        Solver nodes:    [      ^]
-        Indep nodes:     [^      ]
-        Of input nodes, solver nodes for fix_initial are: {1}
-        The indep node is just the first one: {0}
-        """
-        p = brach(transcription='gauss-lobatto', num_segments=1, transcription_order=3,
-                  compressed=True, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0})
-
-    def test_5_lgl(self):
-        """
-        Test one 5th order LGL segment indices
-
-        All nodes are:   [0, 1, 2, 3, 4]
-        Input nodes:     [^     ^     ^]
-        Solver nodes:    [      ^     ^]
-        Of input nodes, solver nodes for fix_initial are: {1, 2}
-        The indep node is just the first one: {0}
-        """
-        p = brach(transcription='gauss-lobatto', num_segments=1, transcription_order=5,
-                  compressed=True, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1, 2})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0})
-
-    def test_3_lgl_compressed(self):
-        """
-        Test one 3rd order LGL segment indices
-
-        All nodes are:   [0, 1, 2, 3, 4]
-        Input nodes:     [^     ^     ^]
-        Solver nodes:    [      ^     ^]
-        Indep nodes:     [^            ]
-        Of input nodes, solver nodes for fix_initial are: {1, 2}
-        The indep node is just the first one: {0}
-        """
-        p = brach(transcription='gauss-lobatto', num_segments=2, transcription_order=3,
-                  compressed=True, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1, 2})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0})
-
-    def test_5_lgl_compressed(self):
-        """
-        Test two 5th order LGL segment indices
-
-        All nodes are:   [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
-        Input nodes:     [^     ^     ^        ^     ^]
-        Solver nodes:    [      ^     ^        ^     ^]
-        Indep nodes:     [^                        ]
-        Of input nodes, solver nodes for fix_initial are: {1, 2, 3, 4}
-        The indep node is just the first one: {0}
-        """
-        p = brach(transcription='gauss-lobatto', num_segments=2, transcription_order=5,
-                  compressed=True, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1, 2, 3, 4})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0})
-
-    def test_3_lgl_uncompressed(self):
-        """
-        Test two 3rd order LGL segment indices
-
-        All nodes are:   [0, 1, 2, 3, 4, 5]
-        Input nodes:     [^     ^  ^     ^]
-        Solver nodes:    [      ^        ^]
-        Indep nodes:     [^        ^      ]
-        Of input nodes, solver nodes for fix_initial are: {1, 3}
-        The indep node is just the first one: {0, 2}
-        """
-        p = brach(transcription='gauss-lobatto', num_segments=2, transcription_order=3,
-                  compressed=False, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1, 3})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0, 2})
-
-    def test_5_lgl_uncompressed(self):
-        """
-        Test two 5th order LGL segment indices
-
-        All nodes are:   [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
-        Input nodes:     [^     ^     ^  ^     ^     ^]
-        Solver nodes:    [      ^     ^        ^     ^]
-        Indep nodes:     [^              ^         ]
-        Of input nodes, solver nodes for fix_initial are: {1, 2, 4, 5}
-        The indep node is just the first one in each segment: {0, 3}
-        """
-        p = brach(transcription='gauss-lobatto', num_segments=2, transcription_order=5,
-                  compressed=False, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1, 2, 4, 5})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0, 3})
-
-    def test_3_radau(self):
-        """
-        Test one 3rd order radau segment indices
-
-        All nodes are:   [0, 1, 2, 3]
-        Input nodes:     [^  ^  ^  ^]
-        Solver nodes:    [   ^  ^  ^]
-        Indep nodes:     [^         ]
-        Of input nodes, solver nodes for fix_initial are: {1, 2, 3}
-        The indep node is just the first one: {0}
-        """
-        p = brach(transcription='radau-ps', num_segments=1, transcription_order=3,
-                  compressed=True, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1, 2, 3})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0})
-
-    def test_5_radau(self):
-        """
-        Test one 5th order radau segment indices
-
-        All nodes are:   [0, 1, 2, 3, 4, 5]
-        Input nodes:     [^  ^  ^  ^  ^  ^]
-        Solver nodes:    [   ^  ^  ^  ^  ^]
-        Indep nodes:     [^         ]
-        Of input nodes, solver nodes for fix_initial are: {1, 2, 3, 4, 5}
-        The indep node is just the first one: {0}
-        """
-        p = brach(transcription='radau-ps', num_segments=1, transcription_order=5,
-                  compressed=True, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1, 2, 3, 4, 5})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0})
-
-    def test_3_radau_compressed(self):
-        """
-        Test one 3rd order radau segment indices
-
-        All nodes are:   [0, 1, 2, 3, 4, 5, 6, 7]
-        Input nodes:     [^  ^  ^  ^     ^  ^  ^]
-        Solver nodes:    [   ^  ^  ^     ^  ^  ^]
-        Indep nodes:     [^            ]
-        Of input nodes, solver nodes for fix_initial are: {1, 2, 3, 4, 5, 6}
-        The indep node is just the first one: {0}
-        """
-        p = brach(transcription='radau-ps', num_segments=2, transcription_order=3,
-                  compressed=True, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1, 2, 3, 4, 5, 6})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0})
-
-    def test_5_radau_compressed(self):
-        """
-        Test two 5th order radau segment indices
-
-        All nodes are:   [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
-        Input nodes:     [^  ^  ^  ^  ^  ^     ^  ^  ^   ^   ^]
-        Solver nodes:    [   ^  ^  ^  ^  ^     ^  ^  ^   ^   ^]
-        Indep nodes:     [^                                   ]
-        Of input nodes, solver nodes for fix_initial are: {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}
-        The indep node is just the first one: {0}
-        """
-        p = brach(transcription='radau-ps', num_segments=2, transcription_order=5,
-                  compressed=True, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1, 2, 3, 4, 5, 6, 7, 8, 9, 10})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0})
-
-    def test_3_radau_uncompressed(self):
-        """
-        Test one 3rd order radau segment indices
-
-        All nodes are:   [0, 1, 2, 3, 4, 5, 6, 7]
-        Input nodes:     [^  ^  ^  ^  ^  ^  ^  ^]
-        Solver nodes:    [   ^  ^  ^     ^  ^  ^]
-        Indep nodes:     [^           ^         ]
-        Of input nodes, solver nodes for fix_initial are: {1, 2, 3, 5, 6, 7}
-        The indep node is just the first one: {0, 4}
-        """
-        p = brach(transcription='radau-ps', num_segments=2, transcription_order=3,
-                  compressed=False, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1, 2, 3, 5, 6, 7})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0, 4})
-
-    def test_5_radau_uncompressed(self):
-        """
-        Test two 5th order radau segment indices
-
-        All nodes are:   [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
-        Input nodes:     [^  ^  ^  ^  ^  ^  ^  ^  ^  ^   ^   ^]
-        Solver nodes:    [   ^  ^  ^  ^  ^     ^  ^  ^   ^   ^]
-        Indep nodes:     [^                 ^                 ]
-        Of input nodes, solver nodes for fix_initial are: {1, 2, 3, 4, 5, 7, 8, 9, 10, 11}
-        The indep node is just the first one: {0, 6}
-        """
-        p = brach(transcription='radau-ps', num_segments=2, transcription_order=5,
-                  compressed=False, solve_segments='forward', run_driver=True)
-
-        state_indeps_comp = p.model.traj0.phases.phase0.indep_states
-
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['solver']), {1, 2, 3, 4, 5, 7, 8, 9, 10, 11})
-        self.assertSetEqual(set(state_indeps_comp.state_idx_map['x']['indep']), {0, 6})
-
-
 @use_tempdirs
 class TestCollocationBalanceApplyNL(unittest.TestCase):
 
@@ -268,9 +28,9 @@ def make_prob(self, transcription, num_segments, transcription_order, compressed
                                 order=transcription_order,
                                 compressed=compressed,)
         elif transcription == 'radau-ps':
-            t = dm.Radau(num_segments=num_segments,
-                         order=transcription_order,
-                         compressed=compressed)
+            t = dm.RadauDeprecated(num_segments=num_segments,
+                                   order=transcription_order,
+                                   compressed=compressed)
 
         traj = dm.Trajectory()
         phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)
diff --git a/dymos/transcriptions/pseudospectral/components/test/test_radau_iter_group.py b/dymos/transcriptions/pseudospectral/components/test/test_radau_iter_group.py
new file mode 100644
index 000000000..4e6158411
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/test/test_radau_iter_group.py
@@ -0,0 +1,229 @@
+import unittest
+
+import numpy as np
+
+import openmdao.api as om
+
+import dymos
+from openmdao.utils.assert_utils import assert_check_partials, assert_near_equal
+from openmdao.utils.testing_utils import use_tempdirs
+
+from dymos.examples.brachistochrone.brachistochrone_vector_states_ode import BrachistochroneVectorStatesODE
+from dymos.utils.misc import GroupWrapperConfig
+from dymos.transcriptions.pseudospectral.components import RadauIterGroup
+from dymos.phase.options import StateOptionsDictionary, TimeOptionsDictionary, ControlOptionsDictionary
+from dymos.transcriptions.grid_data import RadauGrid
+from dymos.utils.testing_utils import PhaseStub, SimpleODE, SimpleVectorizedODE
+
+
+RadauIterGroup = GroupWrapperConfig(RadauIterGroup, [PhaseStub()])
+
+
+@use_tempdirs
+class TestRadauIterGroup(unittest.TestCase):
+
+    def test_solve_segments(self):
+        with dymos.options.temporary(include_check_partials=True):
+            for direction in ['forward', 'backward']:
+                for compressed in [True, False]:
+                    with self.subTest(msg=f'{direction=} {compressed=}'):
+
+                        state_options = {'x': StateOptionsDictionary()}
+
+                        state_options['x']['shape'] = (1,)
+                        state_options['x']['units'] = 's**2'
+                        state_options['x']['targets'] = ['x']
+                        state_options['x']['initial_bounds'] = (None, None)
+                        state_options['x']['final_bounds'] = (None, None)
+                        state_options['x']['solve_segments'] = direction
+                        state_options['x']['rate_source'] = 'x_dot'
+
+                        time_options = TimeOptionsDictionary()
+                        grid_data = RadauGrid(num_segments=2, nodes_per_seg=8, compressed=compressed)
+                        nn = grid_data.subset_num_nodes['all']
+                        ode_class = SimpleODE
+
+                        p = om.Problem()
+                        p.model.add_subsystem('radau', RadauIterGroup(state_options=state_options,
+                                                                      time_options=time_options,
+                                                                      grid_data=grid_data,
+                                                                      ode_class=ode_class))
+
+                        radau = p.model._get_subsystem('radau')
+
+                        radau.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, iprint=-1)
+                        radau.linear_solver = om.DirectSolver()
+
+                        p.setup(force_alloc_complex=True)
+
+                        # Instead of using the TimeComp just transform the node segment taus onto [0, 2]
+                        times = grid_data.node_ptau + 1
+
+                        solution = np.reshape(times**2 + 2 * times + 1 - 0.5 * np.exp(times), (nn, 1))
+
+                        # Each segment is of the same length, so dt_dstau is constant.
+                        # dt_dstau is (tf - t0) / 2.0 / num_seg
+                        p.set_val('radau.dt_dstau', (times[-1] / 2.0 / grid_data.num_segments))
+
+                        if direction == 'forward':
+                            p.set_val('radau.initial_states:x', 0.5)
+                        else:
+                            p.set_val('radau.final_states:x', solution[-1])
+
+                        p.set_val('radau.states:x', 0.0)
+                        p.set_val('radau.ode_all.t', times)
+                        p.set_val('radau.ode_all.p', 1.0)
+
+                        p.run_model()
+
+                        x = p.get_val('radau.states:x')
+                        x_0 = p.get_val('radau.initial_states:x')
+                        x_f = p.get_val('radau.final_states:x')
+
+                        idxs = grid_data.subset_node_indices['state_input']
+                        assert_near_equal(solution[idxs], x, tolerance=1.0E-5)
+                        assert_near_equal(solution[np.newaxis, 0], x_0, tolerance=1.0E-7)
+                        assert_near_equal(solution[np.newaxis, -1], x_f, tolerance=1.0E-7)
+
+                        cpd = p.check_partials(method='cs', compact_print=False, out_stream=None)
+                        assert_check_partials(cpd)
+
+    def test_solve_segments_vector_states(self):
+        with dymos.options.temporary(include_check_partials=True):
+            for direction in ['forward', 'backward']:
+                    with self.subTest(msg=f'{direction=}'):
+
+                        state_options = {'z': StateOptionsDictionary()}
+
+                        state_options['z']['shape'] = (2,)
+                        state_options['z']['units'] = 's**2'
+                        state_options['z']['targets'] = ['z']
+                        state_options['z']['initial_bounds'] = (None, None)
+                        state_options['z']['final_bounds'] = (None, None)
+                        state_options['z']['solve_segments'] = direction
+                        state_options['z']['rate_source'] = 'z_dot'
+
+                        time_options = TimeOptionsDictionary()
+                        grid_data = RadauGrid(num_segments=5, nodes_per_seg=4, compressed=False)
+                        nn = grid_data.subset_num_nodes['all']
+                        ode_class = SimpleVectorizedODE
+
+                        p = om.Problem()
+                        p.model.add_subsystem('radau', RadauIterGroup(state_options=state_options,
+                                                                      time_options=time_options,
+                                                                      grid_data=grid_data,
+                                                                      ode_class=ode_class))
+
+                        radau = p.model._get_subsystem('radau')
+
+                        radau.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, maxiter=100, iprint=-1, atol=1.0E-10, rtol=1.0E-10)
+                        radau.nonlinear_solver.linesearch = om.ArmijoGoldsteinLS(bound_enforcement='vector')
+                        radau.linear_solver = om.DirectSolver()
+
+                        p.setup(force_alloc_complex=True)
+
+                        # Instead of using the TimeComp just transform the node segment taus onto [0, 2]
+                        times = grid_data.node_ptau + 1
+
+                        solution = np.reshape(times**2 + 2 * times + 1 - 0.5 * np.exp(times), (nn, 1))
+
+                        # Each segment is of the same length, so dt_dstau is constant.
+                        # dt_dstau is (tf - t0) / 2.0 / num_seg
+                        p.set_val('radau.dt_dstau', (times[-1] / 2.0 / grid_data.num_segments))
+
+                        if direction == 'forward':
+                            p.set_val('radau.initial_states:z', [0.5, 0.0])
+                        else:
+                            p.set_val('radau.final_states:z', solution[-1], indices=om.slicer[:, 0])
+                            p.set_val('radau.final_states:z', 20.0, indices=om.slicer[:, 1])
+
+                        p.set_val('radau.states:z', 0.0)
+                        p.set_val('radau.ode_all.t', times)
+                        p.set_val('radau.ode_all.p', 1.0)
+
+                        p.run_model()
+
+                        x = p.get_val('radau.states:z')
+                        x_0 = p.get_val('radau.initial_states:z')
+                        x_f = p.get_val('radau.final_states:z')
+
+                        idxs = grid_data.subset_node_indices['state_input']
+                        assert_near_equal(solution[idxs], x[np.newaxis, :, 0].T, tolerance=1.0E-5)
+                        assert_near_equal(solution[0], x_0[:, 0], tolerance=1.0E-5)
+                        assert_near_equal(solution[-1], x_f[:, 0], tolerance=1.0E-5)
+
+                        with np.printoptions(linewidth=1024, edgeitems=1024):
+                            cpd = p.check_partials(method='cs', compact_print=False, show_only_incorrect=True, out_stream=None)
+                        assert_check_partials(cpd, atol=1.0E-5, rtol=1.0)
+
+    def test_autoivc_no_solve(self):
+        with dymos.options.temporary(include_check_partials=True):
+            for direction in ['forward', 'backward']:
+                for compressed in [True, False]:
+                    with self.subTest(msg=f'{direction=} {compressed=}'):
+
+                        state_options = {'x': StateOptionsDictionary()}
+
+                        state_options['x']['shape'] = (1,)
+                        state_options['x']['units'] = 's**2'
+                        state_options['x']['targets'] = ['x']
+                        state_options['x']['initial_bounds'] = (None, None)
+                        state_options['x']['final_bounds'] = (None, None)
+                        state_options['x']['solve_segments'] = False
+                        state_options['x']['rate_source'] = 'x_dot'
+
+                        time_options = TimeOptionsDictionary()
+                        grid_data = RadauGrid(num_segments=2, nodes_per_seg=8, compressed=compressed)
+                        nn = grid_data.subset_num_nodes['all']
+                        ode_class = SimpleODE
+
+                        p = om.Problem()
+                        p.model.add_subsystem('radau', RadauIterGroup(state_options=state_options,
+                                                                      time_options=time_options,
+                                                                      grid_data=grid_data,
+                                                                      ode_class=ode_class))
+
+                        radau = p.model._get_subsystem('radau')
+
+                        p.setup(force_alloc_complex=True)
+
+                        # Instead of using the TimeComp just transform the node segment taus onto [0, 2]
+                        times = grid_data.node_ptau + 1
+
+                        solution = np.reshape(times**2 + 2 * times + 1 - 0.5 * np.exp(times), (nn, 1))
+                        seg_end_idxs = grid_data.subset_node_indices['segment_ends'][1::2][:-1]
+                        solution_input_nodes = np.delete(solution, seg_end_idxs)
+
+                        # Each segment is of the same length, so dt_dstau is constant.
+                        # dt_dstau is (tf - t0) / 2.0 / num_seg
+                        p.set_val('radau.dt_dstau', (times[-1] / 2.0 / grid_data.num_segments))
+
+                        p.set_val('radau.initial_states:x', 0.5)
+                        p.set_val('radau.final_states:x', solution[-1])
+
+                        if compressed:
+                            # Solution is definted at all nodes, need to only provide the input node values.
+                            p.set_val('radau.states:x', solution_input_nodes)
+                        else:
+                            p.set_val('radau.states:x', solution)
+
+                        p.set_val('radau.ode_all.t', times)
+                        p.set_val('radau.ode_all.p', 1.0)
+
+                        p.run_model()
+
+                        x = p.get_val('radau.states:x')
+                        x_0 = p.get_val('radau.initial_states:x')
+                        x_f = p.get_val('radau.final_states:x')
+
+                        idxs = grid_data.subset_node_indices['state_input']
+                        assert_near_equal(solution[idxs], x, tolerance=1.0E-5)
+                        assert_near_equal(solution[np.newaxis, 0], x_0, tolerance=1.0E-7)
+                        assert_near_equal(solution[np.newaxis, -1], x_f, tolerance=1.0E-7)
+
+                        cpd = p.check_partials(method='cs', compact_print=False, out_stream=None)
+                        assert_check_partials(cpd)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/dymos/transcriptions/pseudospectral/gauss_lobatto.py b/dymos/transcriptions/pseudospectral/gauss_lobatto.py
index 6389a50f6..6bd986461 100644
--- a/dymos/transcriptions/pseudospectral/gauss_lobatto.py
+++ b/dymos/transcriptions/pseudospectral/gauss_lobatto.py
@@ -82,7 +82,8 @@ def configure_time(self, phase):
                               src_indices=disc_src_idxs, flat_src_indices=True)
 
         for name, targets in [('t_initial', options['t_initial_targets']),
-                              ('t_duration', options['t_duration_targets'])]:
+                              ('t_duration', options['t_duration_targets']),
+                              ('t_final', options['t_final_targets'])]:
             for t in targets:
                 shape = ode_inputs[t]['shape']
                 if shape == (1,):
diff --git a/dymos/transcriptions/pseudospectral/radau_new.py b/dymos/transcriptions/pseudospectral/radau_new.py
new file mode 100644
index 000000000..5405be9a8
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/radau_new.py
@@ -0,0 +1,934 @@
+import numpy as np
+
+import openmdao.api as om
+
+from ..transcription_base import TranscriptionBase
+from ..common import TimeComp, TimeseriesOutputGroup, TimeseriesOutputComp, ControlComp
+from .components import RadauIterGroup, RadauBoundaryGroup
+
+from ..grid_data import RadauGrid
+from dymos.utils.misc import get_rate_units
+from dymos.utils.introspection import get_promoted_vars, get_source_metadata
+from dymos.utils.indexing import get_constraint_flat_idxs, get_src_indices_by_row
+
+
+class RadauNew(TranscriptionBase):
+    """
+    Radau Pseudospectral Transcription.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+
+    """
+    def __init__(self, **kwargs):
+        super(RadauNew, self).__init__(**kwargs)
+        self._rhs_source = 'ode_iter_group.ode_all'
+
+    def initialize(self):
+        """
+        Declare transcription options.
+        """
+        self.options.declare('grid', types=(RadauGrid, str),
+                             allow_none=True, default=None,
+                             desc='The grid distribution used to layout the control input and interpolation nodes.')
+
+        self.options.declare(name='solve_segments', default=False,
+                             values=(False, 'forward', 'backward'),
+                             desc='Applies \'solve_segments\' behavior to _all_ states in the Phase. '
+                                  'If \'forward\', collocation defects within each '
+                                  'segment are solved with a Newton solver by fixing the initial value in the '
+                                  'phase (if using compressed transcription) or segment (if not using '
+                                  'compressed transcription). This provides a forward shooting (or multiple shooting) '
+                                  'method.  If \'backward\', the final value in the phase or segment is fixed '
+                                  'and a solver finds the other ones to mimic reverse propagation. Set '
+                                  'to False (the default) to explicitly disable the use of a solver to '
+                                  'converge the state time history.')
+
+    def init_grid(self):
+        """
+        Set up the GridData object for the Transcription.
+        """
+        self.grid_data = RadauGrid(num_segments=self.options['num_segments'],
+                                   nodes_per_seg=self.options['order'] + 1,
+                                   segment_ends=self.options['segment_ends'],
+                                   compressed=self.options['compressed'])
+
+    def setup_time(self, phase):
+        """
+        Set up the time component.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        grid_data = self.grid_data
+
+        super(RadauNew, self).setup_time(phase)
+
+        time_comp = TimeComp(num_nodes=grid_data.num_nodes, node_ptau=grid_data.node_ptau,
+                             node_dptau_dstau=grid_data.node_dptau_dstau,
+                             units=phase.time_options['units'],
+                             initial_val=phase.time_options['initial_val'],
+                             duration_val=phase.time_options['duration_val'])
+
+        phase.add_subsystem('time', time_comp, promotes_inputs=['*'], promotes_outputs=['*'])
+
+    def configure_time(self, phase):
+        """
+        Configure the inputs/outputs on the time component.
+
+        This method assumes that target introspection has already been performed by the phase and thus
+        options['targets'], options['time_phase_targets'], options['t_initial_targets'],
+        and options['t_duration_targets'] are all correctly populated.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        super(RadauNew, self).configure_time(phase)
+        phase.time.configure_io()
+        options = phase.time_options
+        ode = phase._get_subsystem(self._rhs_source)
+        ode_inputs = get_promoted_vars(ode, iotypes='input')
+
+        # The tuples here are (name, user_specified_targets, dynamic)
+        for name, targets in [('t', options['targets']),
+                              ('t_phase', options['time_phase_targets'])]:
+            if targets:
+                src_idxs = self.grid_data.subset_node_indices['all']
+                phase.connect(name, [f'ode_all.{t}' for t in targets], src_indices=src_idxs,
+                              flat_src_indices=True)
+                src_idxs = om.slicer[[0, -1], ...]
+                phase.connect(name, [f'boundary_vals.{t}' for t in targets], src_indices=src_idxs)
+
+        for name, targets in [('t_initial', options['t_initial_targets']),
+                              ('t_duration', options['t_duration_targets']),
+                              ('t_final', options['t_final_targets'])]:
+            for t in targets:
+                shape = ode_inputs[t]['shape']
+
+                if shape == (1,):
+                    src_idxs = endpoint_src_idxs = None
+                    flat_src_idxs = None
+                    src_shape = None
+                else:
+                    src_idxs = np.zeros(self.grid_data.subset_num_nodes['all'], dtype=int)
+                    flat_src_idxs = True
+                    src_shape = (1,)
+
+                phase.ode_iter_group.promotes('ode_all', any=[(t, name)], src_indices=src_idxs,
+                                              flat_src_indices=flat_src_idxs, src_shape=src_shape)
+
+                phase.promotes('boundary_vals', inputs=[(t, name)], src_indices=endpoint_src_idxs,
+                               flat_src_indices=flat_src_idxs, src_shape=src_shape)
+            if targets:
+                phase.set_input_defaults(name=name,
+                                         val=np.ones((1,)),
+                                         units=options['units'])
+
+    def setup_states(self, phase):
+        """
+        Set up the states for this transcription.
+
+        In the Radau2 transcription, everything typically done in this
+        method is instead done by the RadauIterGroup.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        self.any_solved_segs = False
+        self.any_connected_opt_segs = False
+        for options in phase.state_options.values():
+            # Transcription solve_segments overrides state solve_segments if its not set
+            if options['solve_segments'] is None:
+                options['solve_segments'] = self.options['solve_segments']
+
+            if options['solve_segments']:
+                self.any_solved_segs = True
+            elif options['input_initial']:
+                self.any_connected_opt_segs = True
+
+    def setup_controls(self, phase):
+        """
+        Setup the control group.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        phase._check_control_options()
+
+        if phase.control_options:
+            gd = self.grid_data
+            control_options = phase.control_options
+            time_units = phase.time_options['units']
+
+            phase.add_subsystem('control_comp',
+                                subsys=ControlComp(time_units=time_units,
+                                                   grid_data=gd,
+                                                   compressed=False,
+                                                   control_options=control_options),
+                                promotes_inputs=['*controls*'],
+                                promotes_outputs=['control_values:*', 'control_rates:*'])
+
+            phase.connect('t_duration_val', 'control_comp.t_duration')
+
+    def configure_controls(self, phase):
+        """
+        Configure the inputs/outputs for the controls.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        super().configure_controls(phase)
+
+        if phase.control_options:
+            phase._get_subsystem('control_comp').configure_io()
+
+            phase.connect('dt_dstau', 'control_comp.dt_dstau')
+
+        for name, options in phase.control_options.items():
+            if options['targets']:
+                phase.connect(f'control_values:{name}', [f'ode_all.{t}' for t in options['targets']])
+                phase.connect(f'control_values:{name}', [f'boundary_vals.{t}' for t in options['targets']],
+                              src_indices=om.slicer[[0, -1], ...])
+
+            if options['rate_targets']:
+                phase.connect(f'control_rates:{name}_rate',
+                              [f'ode_all.{t}' for t in options['rate_targets']])
+                phase.connect(f'control_rates:{name}_rate',
+                              [f'boundary_vals.{t}' for t in options['rate_targets']],
+                              src_indices=om.slicer[[0, -1], ...])
+
+            if options['rate2_targets']:
+                phase.connect(f'control_rates:{name}_rate2',
+                              [f'ode_all.{t}' for t in options['rate2_targets']])
+                phase.connect(f'control_rates:{name}_rate2',
+                              [f'boundary_vals.{t}' for t in options['rate2_targets']],
+                              src_indices=om.slicer[[0, -1], ...])
+
+    def setup_ode(self, phase):
+        """
+        Set up the ode for this transcription.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        ODEClass = phase.options['ode_class']
+        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=RadauIterGroup(grid_data=grid_data, state_options=phase.state_options,
+                                                  time_options=phase.time_options,
+                                                  ode_class=ODEClass,
+                                                  ode_init_kwargs=ode_init_kwargs),
+                            promotes=['*'])
+
+        phase.add_subsystem('boundary_vals',
+                            subsys=RadauBoundaryGroup(grid_data=grid_data,
+                                                      ode_class=ODEClass,
+                                                      ode_init_kwargs=ode_init_kwargs,
+                                                      initial_boundary_constraints=ibcs,
+                                                      final_boundary_constraints=fbcs,
+                                                      objectives=phase._objectives),
+                            promotes_inputs=['initial_states:*', 'final_states:*'])
+
+    def configure_ode(self, phase):
+        """
+        Create connections to the introspected states.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        phase._get_subsystem('boundary_vals').configure_io(phase)
+        phase._get_subsystem('ode_iter_group').configure_io(phase)
+
+    def setup_defects(self, phase):
+        """
+        Pass setup_defects for this transcription.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+        # if any(self._requires_continuity_constraints(phase)):
+        #     phase.add_subsystem('continuity_comp',
+        #                         RadauPSContinuityComp(grid_data=self.grid_data,
+        #                                               # State continuity defects are in the RadauDefectComp.
+        #                                               state_options={},
+        #                                               control_options=phase.control_options,
+        #                                               time_units=phase.time_options['units']))
+
+    def configure_defects(self, phase):
+        """
+        Configure the continuity_comp and connect the collocation constraints.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        grid_data = self.grid_data
+
+        for name, options in phase.state_options.items():
+            rate_source_type = phase.classify_var(options['rate_source'])
+            rate_src_path = self._get_rate_source_path(name, phase)
+            if rate_source_type not in ('state', 'ode'):
+                phase.connect(rate_src_path, f'f_ode:{name}',
+                              src_indices=grid_data.subset_node_indices['col'])
+
+        # any_state_cnty, any_control_cnty, any_control_rate_cnty = self._requires_continuity_constraints(phase)
+
+        # if not any((any_state_cnty, any_control_cnty, any_control_rate_cnty)):
+        #     return
+
+        # # Add the continuity constraint component if necessary
+        # segment_end_idxs = grid_data.subset_node_indices['segment_ends']
+
+        # if any_control_rate_cnty:
+        #     phase.connect('t_duration_val', 'continuity_comp.t_duration')
+
+        # phase.continuity_comp.configure_io()
+
+        # for name, options in phase.control_options.items():
+        #     # The sub-indices of control_disc indices that are segment ends
+        #     segment_end_idxs = grid_data.subset_node_indices['segment_ends']
+        #     src_idxs = get_src_indices_by_row(segment_end_idxs, options['shape'], flat=True)
+
+        #     # enclose indices in tuple to ensure shaping of indices works
+        #     src_idxs = (src_idxs,)
+
+        #     if options['continuity']:
+        #         phase.connect(f'control_values:{name}',
+        #                       f'continuity_comp.controls:{name}',
+        #                       src_indices=src_idxs, flat_src_indices=True)
+
+        #     if options['rate_continuity']:
+        #         phase.connect(f'control_rates:{name}_rate',
+        #                       f'continuity_comp.control_rates:{name}_rate',
+        #                       src_indices=src_idxs, flat_src_indices=True)
+
+        #     if options['rate2_continuity']:
+        #         phase.connect(f'control_rates:{name}_rate2',
+        #                       f'continuity_comp.control_rates:{name}_rate2',
+        #                       src_indices=src_idxs, flat_src_indices=True)
+
+    def setup_duration_balance(self, phase):
+        """
+        Set up the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
+    def configure_duration_balance(self, phase):
+        """
+        Configure the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
+    def setup_solvers(self, phase):
+        """
+        Set up the solvers.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
+    def configure_solvers(self, phase, requires_solvers=None):
+        """
+        Configure the solvers.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        requires_solvers : dict[str: bool]
+            A dictionary mapping a string descriptor of a reason why a solver is required,
+            and whether a solver is required.
+        """
+        ode_iter_group = phase._get_subsystem('ode_iter_group')
+        req_solvers = {'implicit outputs': ode_iter_group._implicit_outputs}
+
+        if requires_solvers is not None:
+            req_solvers.update(requires_solvers)
+
+        super().configure_solvers(phase, requires_solvers=req_solvers)
+
+    def configure_timeseries_outputs(self, phase):
+        """
+        Create connections from time series to all post-introspection sources.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        gd = self.grid_data
+        nsin = gd.subset_num_nodes['state_input']
+        state_options = phase.state_options
+        for timeseries_name, timeseries_options in phase._timeseries.items():
+            timeseries_comp = phase._get_subsystem(f'{timeseries_name}.timeseries_comp')
+            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,
+                # or connect it otherwise.
+                if src.startswith('states:'):
+                    state_name = src.split(':')[-1]
+                    ts_inputs_to_promote.append((input_name, f'states:{state_name}'))
+                    src_shape = (nsin,) + state_options[state_name]['shape']
+                    phase.promotes(timeseries_name,
+                                   inputs=[(input_name, f'states:{state_name}')],
+                                   src_shape=src_shape, src_indices=src_idxs)
+                else:
+                    phase.connect(src_name=src,
+                                  tgt_name=f'{timeseries_name}.{input_name}',
+                                  src_indices=src_idxs)
+
+    def setup_timeseries_outputs(self, phase):
+        """
+        Set up the timeseries for this transcription.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        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:
+                ogd = options['transcription'].grid_data
+
+            timeseries_comp = TimeseriesOutputComp(input_grid_data=gd,
+                                                   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.connect('dt_dstau', f'{name}.dt_dstau', flat_src_indices=True)
+
+    def _get_constraint_kwargs(self, constraint_type, options, phase):
+        """
+        Given the constraint options provide the keyword arguments for the OpenMDAO add_constraint method.
+
+        Parameters
+        ----------
+        constraint_type : str
+            One of 'initial', 'final', or 'path'.
+        options : dict
+            The constraint options.
+        phase : Phase
+            The dymos phase to which the constraint applies.
+
+        Returns
+        -------
+        con_output : str
+            The phase-relative path being constrained.
+        constraint_kwargs : dict
+            Keyword arguments for the OpenMDAO add_constraint method.
+        """
+        num_nodes = self._get_num_timeseries_nodes()
+
+        constraint_kwargs = {key: options for key, options in options.items()}
+        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_type = phase.classify_var(var)
+
+        # These are the flat indices at a single point in time used
+        # in either initial, final, or path constraints.
+        idxs_in_initial = phase._indices_in_constraints(var, 'initial')
+        idxs_in_final = phase._indices_in_constraints(var, 'final')
+        idxs_in_path = phase._indices_in_constraints(var, 'path')
+
+        size = np.prod(options['shape'], dtype=int)
+
+        flat_idxs = get_constraint_flat_idxs(options)
+
+        # Now we need to convert the indices given by the user at any given point
+        # to flat indices to be given to OpenMDAO as flat indices spanning the phase.
+        if var_type == 'parameter':
+            if any([idxs_in_initial.intersection(idxs_in_final),
+                    idxs_in_initial.intersection(idxs_in_path),
+                    idxs_in_final.intersection(idxs_in_path)]):
+                raise RuntimeError(f'In phase {phase.pathname}, parameter `{var}` is subject to multiple boundary '
+                                   f'or path constraints.\nParameters are single values that do not change in '
+                                   f'time, and may only be used in a single boundary or path constraint.')
+            constraint_kwargs['indices'] = flat_idxs
+        else:
+            if constraint_type == 'initial':
+                constraint_kwargs['indices'] = flat_idxs
+            elif constraint_type == 'final':
+                constraint_kwargs['indices'] = size + flat_idxs
+            else:
+                # Path
+                path_idxs = []
+                for i in range(num_nodes):
+                    path_idxs.extend(size * i + flat_idxs)
+
+                constraint_kwargs['indices'] = path_idxs
+
+        alias_map = {'path': 'path_constraint',
+                     'initial': 'initial_boundary_constraint',
+                     'final': 'final_boundary_constraint'}
+
+        str_idxs = '' if options['indices'] is None else f'{options["indices"]}'
+
+        constraint_kwargs['alias'] = f'{phase.pathname}->{alias_map[constraint_type]}->{con_name}{str_idxs}'
+        constraint_kwargs.pop('name')
+        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
+
+    def _get_objective_src(self, var, loc, phase, ode_outputs=None):
+        """
+        Return the path to the variable that will be used as the objective.
+
+        Parameters
+        ----------
+        var : str
+            Name of the variable to be used as the objective.
+        loc : str
+            The location of the objective in the phase ['initial', 'final'].
+        phase : dymos.Phase
+            Phase object containing in which the objective resides.
+        ode_outputs : dict or None
+            A dictionary of ODE outputs as returned by get_promoted_vars.
+
+        Returns
+        -------
+        obj_path : str
+            Path to the source.
+        shape : tuple
+            Source shape.
+        units : str
+            Source units.
+        linear : bool
+            True if the objective quantity1 is linear.
+        """
+        time_units = phase.time_options['units']
+        var_type = phase.classify_var(var)
+
+        if ode_outputs is None:
+            ode_outputs = get_promoted_vars(phase._get_subsystem(self._rhs_source), 'output')
+
+        if var_type == 't':
+            shape = (1,)
+            units = time_units
+            linear = True
+            constraint_path = 't'
+        elif var_type == 't_phase':
+            shape = (1,)
+            units = time_units
+            linear = True
+            constraint_path = 't_phase'
+        elif var_type == 'state':
+            shape = phase.state_options[var]['shape']
+            units = phase.state_options[var]['units']
+            linear = False
+            constraint_path = f'boundary_vals.{var}'
+        elif var_type == 'indep_control':
+            shape = phase.control_options[var]['shape']
+            units = phase.control_options[var]['units']
+            linear = True
+            constraint_path = f'control_values:{var}'
+        elif var_type == 'input_control':
+            shape = phase.control_options[var]['shape']
+            units = phase.control_options[var]['units']
+            linear = False
+            constraint_path = f'control_values:{var}'
+        elif var_type == 'indep_polynomial_control':
+            shape = phase.control_options[var]['shape']
+            units = phase.control_options[var]['units']
+            linear = True
+            constraint_path = f'control_values:{var}'
+        elif var_type == 'input_polynomial_control':
+            shape = phase.control_options[var]['shape']
+            units = phase.control_options[var]['units']
+            linear = False
+            constraint_path = f'control_values:{var}'
+        elif var_type == 'parameter':
+            shape = phase.parameter_options[var]['shape']
+            units = phase.parameter_options[var]['units']
+            linear = True
+            constraint_path = f'parameter_vals:{var}'
+        elif var_type in ('control_rate', 'control_rate2'):
+            control_var = var[:-5] if var_type == 'control_rate' else var[:-6]
+            shape = phase.control_options[control_var]['shape']
+            control_units = phase.control_options[control_var]['units']
+            d = 2 if var_type == 'control_rate2' else 1
+            control_rate_units = get_rate_units(control_units, time_units, deriv=d)
+            units = control_rate_units
+            linear = False
+            constraint_path = f'control_rates:{var}'
+        elif var_type in ('polynomial_control_rate', 'polynomial_control_rate2'):
+            control_var = var[:-5]
+            shape = phase.control_options[control_var]['shape']
+            control_units = phase.control_options[control_var]['units']
+            d = 2 if var_type == 'polynomial_control_rate2' else 1
+            control_rate_units = get_rate_units(control_units, time_units, deriv=d)
+            units = control_rate_units
+            linear = False
+            constraint_path = f'control_rates:{var}'
+        elif var_type == 'timeseries_exec_comp_output':
+            shape = (1,)
+            units = None
+            constraint_path = f'timeseries.timeseries_exec_comp.{var}'
+            linear = False
+        else:
+            # Failed to find variable, assume it is in the ODE. This requires introspection.
+            constraint_path = f'boundary_vals.{var}'
+            meta = get_source_metadata(ode_outputs, var, user_units=None, user_shape=None)
+            shape = meta['shape']
+            units = meta['units']
+            linear = False
+
+        return constraint_path, shape, units, linear
+
+    def _get_num_timeseries_nodes(self):
+        """
+        Return the number of nodes in the default timeseries for this transcription.
+
+        Returns
+        -------
+        int
+            The number of nodes in the default timeseries for this transcription.
+        """
+        return self.grid_data.num_nodes
+
+    def _get_rate_source_path(self, state_name, phase):
+        """
+        Return the rate source location and indices for a given state name.
+
+        Parameters
+        ----------
+        state_name : str
+            Name of the state.
+        phase : dymos.Phase
+            Phase object containing the rate source.
+
+        Returns
+        -------
+        str
+            Path to the rate source.
+        ndarray
+            Array of source indices.
+        """
+        try:
+            var = phase.state_options[state_name]['rate_source']
+        except RuntimeError:
+            raise ValueError(f"state '{state_name}' in phase '{phase.name}' was not given a "
+                             "rate_source")
+
+        # Note the rate source must be shape-compatible with the state
+        var_type = phase.classify_var(var)
+
+        # Determine the path to the variable
+        if var_type == 't':
+            rate_path = 't'
+        elif var_type == 't_phase':
+            rate_path = 't_phase'
+        elif var_type == 'state':
+            rate_path = f'states:{var}'
+        elif var_type == 'indep_control':
+            rate_path = f'control_values:{var}'
+        elif var_type == 'input_control':
+            rate_path = f'control_values:{var}'
+        elif var_type == 'control_rate':
+            control_name = var[:-5]
+            rate_path = f'control_rates:{control_name}_rate'
+        elif var_type == 'control_rate2':
+            control_name = var[:-6]
+            rate_path = f'control_rates:{control_name}_rate2'
+        elif var_type == 'indep_polynomial_control':
+            rate_path = f'control_values:{var}'
+        elif var_type == 'input_polynomial_control':
+            rate_path = f'control_values:{var}'
+        elif var_type == 'polynomial_control_rate':
+            control_name = var[:-5]
+            rate_path = f'control_rates:{control_name}_rate'
+        elif var_type == 'polynomial_control_rate2':
+            control_name = var[:-6]
+            rate_path = f'control_rates:{control_name}_rate2'
+        elif var_type == 'parameter':
+            rate_path = f'parameter_vals:{var}'
+        else:
+            # Failed to find variable, assume it is in the ODE
+            rate_path = f'ode_all.{var}'
+
+        return rate_path
+
+    def _get_timeseries_var_source(self, var, output_name, phase):
+        """
+        Return the source path and indices for a given variable to be connected to a timeseries.
+
+        Parameters
+        ----------
+        var : str
+            Name of the timeseries variable whose source is desired.
+        output_name : str
+            Name of the timeseries output whose source is desired.
+        phase : dymos.Phase
+            Phase object containing the variable, either as state, time, control, etc., or as an ODE output.
+
+        Returns
+        -------
+        meta : dict
+            Metadata pertaining to the variable at the given path. This dict contains 'src' (the path to the
+            timeseries source), 'src_idxs' (an array of the
+            source indices), 'units' (the units of the source variable), and 'shape' (the shape of the variable at
+            a given node).
+        """
+        gd = self.grid_data
+        var_type = phase.classify_var(var)
+        time_units = phase.time_options['units']
+
+        transcription = phase.options['transcription']
+        ode = transcription._get_ode(phase)
+        ode_outputs = get_promoted_vars(ode, 'output')
+
+        # The default for node_idxs, applies to everything except states and parameters.
+        node_idxs = gd.subset_node_indices['all']
+
+        meta = {}
+
+        # Determine the path to the variable
+        if var_type == 't':
+            path = 't'
+            src_units = time_units
+            src_shape = (1,)
+        elif var_type == 't_phase':
+            path = 't_phase'
+            src_units = time_units
+            src_shape = (1,)
+        elif var_type == 'state':
+            path = f'states:{var}'
+            src_units = phase.state_options[var]['units']
+            src_shape = phase.state_options[var]['shape']
+
+            # Find the state_input indices which occur at segment endpoints, and repeat them twice
+            state_input_idxs = gd.subset_node_indices['state_input']
+            repeat_idxs = np.ones_like(state_input_idxs)
+            if self.options['compressed']:
+                segment_end_idxs = gd.subset_node_indices['segment_ends'][1:-1]
+                # Repeat nodes that are on segment bounds (but not the first or last nodes in the phase)
+                nodes_to_repeat = list(set(state_input_idxs).intersection(set(segment_end_idxs)))
+                # Now find these nodes in the state input indices
+                idxs_of_ntr_in_state_inputs = np.where(np.isin(state_input_idxs, nodes_to_repeat))[0]
+                # All state input nodes are used once, but nodes_to_repeat are used twice
+                repeat_idxs[idxs_of_ntr_in_state_inputs] = 2
+            # Now we have a way of mapping the state input indices to all nodes
+            map_input_node_idxs_to_all = np.repeat(np.arange(gd.subset_num_nodes['state_input'],
+                                                   dtype=int), repeats=repeat_idxs)
+            # Now select the subset of nodes we want to use.
+            node_idxs = map_input_node_idxs_to_all[gd.subset_node_indices['all']]
+        elif var_type in ['indep_control', 'input_control']:
+            path = f'control_values:{var}'
+            src_units = phase.control_options[var]['units']
+            src_shape = phase.control_options[var]['shape']
+        elif var_type == 'control_rate':
+            control_name = var[:-5]
+            path = f'control_rates:{control_name}_rate'
+            control_name = var[:-5]
+            src_units = get_rate_units(phase.control_options[control_name]['units'], time_units, deriv=1)
+            src_shape = phase.control_options[control_name]['shape']
+        elif var_type == 'control_rate2':
+            control_name = var[:-6]
+            path = f'control_rates:{control_name}_rate2'
+            src_units = get_rate_units(phase.control_options[control_name]['units'], time_units, deriv=2)
+            src_shape = phase.control_options[control_name]['shape']
+        elif var_type in ['indep_polynomial_control', 'input_polynomial_control']:
+            path = f'control_values:{var}'
+            src_units = phase.control_options[var]['units']
+            src_shape = phase.control_options[var]['shape']
+        elif var_type == 'polynomial_control_rate':
+            control_name = var[:-5]
+            path = f'control_rates:{control_name}_rate'
+            control = phase.control_options[control_name]
+            src_units = get_rate_units(control['units'], time_units, deriv=1)
+            src_shape = control['shape']
+        elif var_type == 'polynomial_control_rate2':
+            control_name = var[:-6]
+            path = f'control_rates:{control_name}_rate2'
+            control = phase.control_options[control_name]
+            src_units = get_rate_units(control['units'], time_units, deriv=2)
+            src_shape = control['shape']
+        elif var_type == 'parameter':
+            path = f'parameter_vals:{var}'
+            # Timeseries are never a static_target
+            node_idxs = np.zeros(gd.subset_num_nodes['all'], dtype=int)
+            src_units = phase.parameter_options[var]['units']
+            src_shape = phase.parameter_options[var]['shape']
+        else:
+            # Failed to find variable, assume it is in the ODE
+            path = f'ode_all.{var}'
+            meta = get_source_metadata(ode_outputs, src=var)
+            src_shape = meta['shape']
+            src_units = meta['units']
+            src_tags = meta['tags']
+            if 'dymos.static_output' in src_tags:
+                raise RuntimeError(
+                    f'ODE output {var} is tagged with "dymos.static_output" and cannot be a timeseries output.')
+
+        src_idxs = om.slicer[node_idxs, ...]
+
+        meta['src'] = path
+        meta['src_idxs'] = src_idxs
+        meta['units'] = src_units
+        meta['shape'] = src_shape
+
+        return meta
+
+    def get_parameter_connections(self, name, phase):
+        """
+        Return info about a parameter's target connections in the phase.
+
+        Parameters
+        ----------
+        name : str
+            Parameter name.
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+
+        Returns
+        -------
+        list of (paths, indices)
+            A list containing a tuple of target paths and corresponding src_indices to which the
+            given design variable is to be connected.
+        """
+        connection_info = []
+        if name in phase.parameter_options:
+            options = phase.parameter_options[name]
+            for tgt in options['targets']:
+                if tgt in options['static_targets']:
+                    src_idxs = np.squeeze(get_src_indices_by_row([0], options['shape']), axis=0)
+                    endpoint_src_idxs = om.slicer[:, ...]
+                else:
+                    src_idxs_raw = np.zeros(self.grid_data.subset_num_nodes['all'], dtype=int)
+                    src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape'])
+                    endpoint_src_idxs = om.slicer[[0, -1], ...]
+                    if options['shape'] == (1,):
+                        src_idxs = src_idxs.ravel()
+
+                connection_info.append((f'ode_all.{tgt}', (src_idxs,)))
+                connection_info.append((f'boundary_vals.{tgt}', (endpoint_src_idxs,)))
+
+        return connection_info
+
+    def _requires_continuity_constraints(self, phase):
+        """
+        Test whether state and/or control and/or control rate continuity are required.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase to which this transcription applies.
+
+        Returns
+        -------
+        any_state_continuity : bool
+            True if any state continuity is required to be enforced.
+        any_control_continuity : bool
+            True if any control value continuity is required to be enforced.
+        any_control_rate_continuity : bool
+            True if any control rate continuity is required to be enforced.
+
+        """
+        num_seg = self.grid_data.num_segments
+        compressed = self.grid_data.compressed
+
+        any_state_continuity = num_seg > 1 and not compressed
+        any_control_continuity = any([opts['continuity'] for opts in phase.control_options.values()])
+        any_control_continuity = any_control_continuity and num_seg > 1
+        any_rate_continuity = any([opts['rate_continuity'] or opts['rate2_continuity']
+                                   for opts in phase.control_options.values()])
+        any_rate_continuity = any_rate_continuity and num_seg > 1
+
+        return any_state_continuity, any_control_continuity, any_rate_continuity
+
+    def _phase_set_state_val(self, phase, name, vals, time_vals, interpolation_kind):
+        """
+        Provide variable names and values when phase.set_state_vals is called.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase to which this transcription applies.
+        name : str
+            The name of the phase variable to be set.
+        vals : ndarray or Sequence or float
+            Array of control/state/parameter values.
+        time_vals : ndarray or Sequence or None
+            Array of integration variable values.
+        interpolation_kind : str
+            Specifies the kind of interpolation, as per the scipy.interpolate package.
+            One of ('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic'
+            where 'zero', 'slinear', 'quadratic' and 'cubic' refer to a spline
+            interpolation of zeroth, first, second or third order) or as an
+            integer specifying the order of the spline interpolator to use.
+            Default is 'linear'.
+
+        Returns
+        -------
+        input_data : dict
+            Dict containing the values that need to be set in the phase
+
+        """
+        input_data = {}
+        if np.isscalar(vals):
+            input_data[f'states:{name}'] = vals
+            input_data[f'initial_states:{name}'] = vals
+            input_data[f'final_states:{name}'] = vals
+        else:
+            interp_vals = phase.interp(name, vals, time_vals,
+                                       nodes='state_input',
+                                       kind=interpolation_kind)
+            input_data[f'states:{name}'] = interp_vals
+            input_data[f'initial_states:{name}'] = vals[0]
+            input_data[f'final_states:{name}'] = vals[-1]
+
+        return input_data
diff --git a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
index 4457440e7..52f7d9850 100644
--- a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
+++ b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
@@ -65,7 +65,8 @@ def configure_time(self, phase):
                               flat_src_indices=True if dynamic else None)
 
         for name, targets in [('t_initial', options['t_initial_targets']),
-                              ('t_duration', options['t_duration_targets'])]:
+                              ('t_duration', options['t_duration_targets']),
+                              ('t_final', options['t_final_targets'])]:
             for t in targets:
                 shape = ode_inputs[t]['shape']
 
diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index cbced1fb5..3a124d47f 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -32,8 +32,9 @@ def __init__(self, **kwargs):
                              allow_none=True, desc='Locations of segment ends or None for equally '
                              'spaced segments')
         self.options.declare('order', default=3, types=(int, Sequence, np.ndarray),
-                             desc='Order of the state transcription. The order of the control '
-                                  'transcription is `order - 1`.')
+                             desc='Order of the state transcription.',
+                             deprecation='order is deprecated. Use nodes_per_seg to specify the number of '
+                             'nodes in each segment.')
         self.options.declare('compressed', default=True, types=bool,
                              desc='Use compressed transcription, meaning state and control values'
                                   'at segment boundaries are not duplicated on input.  This '
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index f34ceb902..4f67de4c7 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -160,7 +160,9 @@ def _configure_constraint_introspection(phase):
         The phase object whose boundary and path constraints are to be introspected.
     """
     from ..transcriptions import Birkhoff
+    from ..transcriptions.pseudospectral.radau_new import RadauNew
     birkhoff = isinstance(phase.options['transcription'], Birkhoff)
+    radau_new = isinstance(phase.options['transcription'], RadauNew)
 
     for constraint_type, constraints in [('initial', phase._initial_boundary_constraints),
                                          ('final', phase._final_boundary_constraints),
@@ -194,7 +196,7 @@ def _configure_constraint_introspection(phase):
                 state_units = phase.state_options[var]['units']
                 con['shape'] = state_shape
                 con['units'] = state_units if con['units'] is None else con['units']
-                if birkhoff and constraint_type in ('initial', 'final'):
+                if (birkhoff or radau_new) and constraint_type in ('initial', 'final'):
                     con['constraint_path'] = f'boundary_vals.{var}'
                 else:
                     con['constraint_path'] = f'timeseries.{prefix}{var}'
@@ -226,7 +228,7 @@ def _configure_constraint_introspection(phase):
                 con['shape'] = control_shape
                 con['units'] = get_rate_units(control_units, time_units, deriv=1) \
                     if con['units'] is None else con['units']
-                if birkhoff and constraint_type in ('initial', 'final'):
+                if (birkhoff or radau_new) and constraint_type in ('initial', 'final'):
                     con['constraint_path'] = f'boundary_vals.{var}'
                 else:
                     con['constraint_path'] = f'timeseries.{prefix}{var}'
@@ -239,7 +241,7 @@ def _configure_constraint_introspection(phase):
                 con['shape'] = control_shape
                 con['units'] = get_rate_units(control_units, time_units, deriv=2) \
                     if con['units'] is None else con['units']
-                if birkhoff and constraint_type in ('initial', 'final'):
+                if (birkhoff or radau_new) and constraint_type in ('initial', 'final'):
                     con['constraint_path'] = f'boundary_vals.{var}'
                 else:
                     con['constraint_path'] = f'timeseries.{prefix}{var}'
@@ -258,7 +260,7 @@ def _configure_constraint_introspection(phase):
                 con['shape'] = meta['shape']
                 con['units'] = meta['units']
 
-                if birkhoff and constraint_type in ('initial', 'final'):
+                if (birkhoff or radau_new) and constraint_type in ('initial', 'final'):
                     con['constraint_path'] = f'boundary_vals.{var}'
                 else:
                     con['constraint_path'] = f'timeseries.{con["constraint_name"]}'
diff --git a/dymos/utils/testing_utils.py b/dymos/utils/testing_utils.py
index 2d41476ad..54f3f5d14 100644
--- a/dymos/utils/testing_utils.py
+++ b/dymos/utils/testing_utils.py
@@ -410,3 +410,82 @@ def wrap(*args, **kwargs):
                         os.environ[k] = v
 
         return wrap
+
+class PhaseStub():
+    """
+    A stand-in for the Phase during config_io for testing.
+    It just supports the classify_var method and returns "ode", the only value needed for unittests.
+    """
+    def __init__(self):
+        self.nonlinear_solver = None
+        self.linear_solver = None
+
+    def classify_var(self, name):
+        return 'ode'
+
+class SimpleODE(om.ExplicitComponent):
+    """
+    A simple ODE from https://math.okstate.edu/people/yqwang/teaching/math4513_fall11/Notes/rungekutta.pdf
+    """
+    def initialize(self):
+        self.options.declare('num_nodes', types=(int,))
+
+    def setup(self):
+        nn = self.options['num_nodes']
+        self.add_input('x', shape=(nn,), units='s**2')
+        self.add_input('t', shape=(nn,), units='s')
+        self.add_input('p', shape=(nn,), units='s**2')
+
+        self.add_output('x_dot', shape=(nn,), units='s')
+
+        ar = np.arange(nn, dtype=int)
+        self.declare_partials(of='x_dot', wrt='x', rows=ar, cols=ar, val=1.0)
+        self.declare_partials(of='x_dot', wrt='t', rows=ar, cols=ar)
+        self.declare_partials(of='x_dot', wrt='p', rows=ar, cols=ar, val=1.0)
+
+    def compute(self, inputs, outputs):
+        x = inputs['x']
+        t = inputs['t']
+        p = inputs['p']
+        outputs['x_dot'] = x - t**2 + p
+
+    def compute_partials(self, inputs, partials):
+        t = inputs['t']
+        partials['x_dot', 't'] = -2*t
+
+
+class SimpleVectorizedODE(om.ExplicitComponent):
+    """
+    A simple ODE from https://math.okstate.edu/people/yqwang/teaching/math4513_fall11/Notes/rungekutta.pdf
+    """
+    def initialize(self):
+        self.options.declare('num_nodes', types=(int,))
+
+    def setup(self):
+        nn = self.options['num_nodes']
+        self.add_input('z', shape=(nn, 2), units='s**2')
+        self.add_input('t', shape=(nn,), units='s')
+        self.add_input('p', shape=(nn,), units='s**2')
+
+        self.add_output('z_dot', shape=(nn, 2), units='s')
+
+        cs = np.repeat(np.arange(nn, dtype=int), 2)
+        ar2 = np.arange(2 * nn, dtype=int)
+        dzdot_dz_pattern = np.arange(2 * nn, step=2, dtype=int)
+        self.declare_partials(of='z_dot', wrt='z', rows=dzdot_dz_pattern, cols=dzdot_dz_pattern, val=1.0)
+        self.declare_partials(of='z_dot', wrt='t', rows=ar2, cols=cs)
+        dzdot_dp_rows = np.arange(2 * nn, step=2, dtype=int)
+        dzdot_dp_cols = np.arange(nn, dtype=int)
+        self.declare_partials(of='z_dot', wrt='p', rows=dzdot_dp_rows, cols=dzdot_dp_cols, val=1.0)
+
+    def compute(self, inputs, outputs):
+        z = inputs['z']
+        t = inputs['t']
+        p = inputs['p']
+        outputs['z_dot'][:, 0] = z[:, 0] - t**2 + p
+        outputs['z_dot'][:, 1] = 10 * t
+
+    def compute_partials(self, inputs, partials):
+        t = inputs['t']
+        partials['z_dot', 't'][0::2] = -2*t
+        partials['z_dot', 't'][1::2] = 10
\ No newline at end of file
diff --git a/joss/test/test_cannonball_for_joss.py b/joss/test/test_cannonball_for_joss.py
index e47014806..c1e30dc1a 100644
--- a/joss/test/test_cannonball_for_joss.py
+++ b/joss/test/test_cannonball_for_joss.py
@@ -136,7 +136,7 @@ def compute(self, inputs, outputs):
         traj.add_parameter('mass', units='kg', val=0.01, static_target=True)
         traj.add_parameter('area', units='m**2', static_target=True)
 
-        tx = dm.Radau(num_segments=5, order=3, compressed=True)
+        tx = dm.Radau(num_segments=5, order=3, compressed=False)
         ascent = dm.Phase(transcription=tx, ode_class=CannonballODE)
         traj.add_phase('ascent', ascent)
 
@@ -154,7 +154,7 @@ def compute(self, inputs, outputs):
         #     we will set it to zero so that the phase ends at apogee)
         ascent.add_state('gam', units='rad', rate_source='gam_dot',
                          fix_initial=False, fix_final=True)
-        ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100),
+        ascent.set_time_options(fix_initial=True, duration_bounds=(10, 100),
                                 duration_ref=100, units='s')
 
         ascent.add_parameter('mass', units='kg', val=0.01, static_target=True)
@@ -204,23 +204,41 @@ def compute(self, inputs, outputs):
 
         # Set Initial guesses for static dvs and ascent
         p.set_val('static_calcs.radius', 0.05, units='m')
-        p.set_val('traj.ascent.t_duration', 10.0)
 
-        p.set_val('traj.ascent.states:r', ascent.interp('r', [0, 100]))
-        p.set_val('traj.ascent.states:h', ascent.interp('h', [0, 100]))
-        p.set_val('traj.ascent.states:v', ascent.interp('v', [200, 150]))
-        p.set_val('traj.ascent.states:gam', ascent.interp('gam', [25, 0]), units='deg')
+        ascent.set_time_val(initial=0, duration=10.0)
+        ascent.set_state_val('r', [0, 100])
+        ascent.set_state_val('h', [0, 100])
+        ascent.set_state_val('v', [200, 150])
+        ascent.set_state_val('gam', [25, 0], units='deg')
 
-        # more intitial guesses for descent
-        p.set_val('traj.descent.t_initial', 10.0)
-        p.set_val('traj.descent.t_duration', 10.0)
+        # These set_val commands applied to previous versions of dymos.
+        # With the updated transcriptions as of dymos 1.14.0, we would also need to set
+        # values for `traj.adscent.initial_states:{r|h|v}` and
+        # `traj.adscent.final_states:{r|h|v}`.  By using the `set_state_val` interface instead,
+        # dymos internally sets all appropriate variables pertaining to the states.
 
-        p.set_val('traj.descent.states:r', descent.interp('r', [100, 200]))
-        p.set_val('traj.descent.states:h', descent.interp('h', [100, 0]))
-        p.set_val('traj.descent.states:v', descent.interp('v', [150, 200]))
-        p.set_val('traj.descent.states:gam', descent.interp('gam', [0, -45]), units='deg')
+        # p.set_val('traj.ascent.t_duration', 10.0)
+        # p.set_val('traj.ascent.states:r', ascent.interp('r', [0, 100]))
+        # p.set_val('traj.ascent.states:h', ascent.interp('h', [0, 100]))
+        # p.set_val('traj.ascent.states:v', ascent.interp('v', [200, 150]))
+        # p.set_val('traj.ascent.states:gam', ascent.interp('gam', [25, 0]), units='deg')
 
-        dm.run_problem(p, simulate=True, make_plots=True)
+        # more intitial guesses for descent
+        descent.set_time_val(initial=10.0, duration=10.0)
+        descent.set_state_val('r', [100, 200])
+        descent.set_state_val('h', [100, 0])
+        descent.set_state_val('v', [150, 200])
+        descent.set_state_val('gam', [0, -45], units='deg')
+
+        # Pre-dymos 1.14.0 code
+        # p.set_val('traj.descent.t_initial', 10.0)
+        # p.set_val('traj.descent.t_duration', 10.0)
+        # p.set_val('traj.descent.states:r', descent.interp('r', [100, 200]))
+        # p.set_val('traj.descent.states:h', descent.interp('h', [100, 0]))
+        # p.set_val('traj.descent.states:v', descent.interp('v', [150, 200]))
+        # p.set_val('traj.descent.states:gam', descent.interp('gam', [0, -45]), units='deg')
+
+        dm.run_problem(p, run_driver=True, simulate=True, make_plots=True)
 
         fig, ax = plt.subplots()
         x0 = p.get_val('traj.ascent.timeseries.r', units='m')
@@ -236,7 +254,7 @@ def compute(self, inputs, outputs):
         fig.savefig('cannonball_hr.png', bbox_inches='tight')
         # End code for paper
 
-        assert_near_equal(x1[-1], 3064, tolerance=1.0E-4)
+        assert_near_equal(x1[-1], 3064, tolerance=1.0E-2)
 
 
 if __name__ == '__main__':