diff --git a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb
index 98da7d170..22d3b3bd2 100644
--- a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb
+++ b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb
@@ -4,6 +4,7 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
+    "scrolled": true,
     "tags": [
      "active-ipynb",
      "remove-input",
@@ -76,6 +77,7 @@
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
+    "scrolled": true,
     "tags": [
      "output_scroll"
     ]
@@ -148,7 +150,9 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "class CannonballODE(om.ExplicitComponent):\n",
@@ -232,7 +236,9 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [],
    "source": [
     "def initial_guess(t_dur, gam0=np.pi/3):\n",
@@ -265,10 +271,12 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": false
+   },
    "outputs": [],
    "source": [
-    "p = om.Problem(model=om.Group())\n",
+    "p = om.Problem(name='cannonball_implicit_duration', model=om.Group())\n",
     "\n",
     "p.driver = om.pyOptSparseDriver()\n",
     "p.driver.options['optimizer'] = 'IPOPT'\n",
@@ -321,6 +329,14 @@
     "# A nonlinear solver is used to find the duration of required to satisfy the condition.\n",
     "phase.set_duration_balance('h', val=0.0)\n",
     "\n",
+    "# Prior to setup, this phase has the default NonlinearRunOnce and LinearRunOnce solvers.\n",
+    "# Because of the implicit behavior introduced by set_duration_balance, it will automatically\n",
+    "# use a NewtonSolver with an Armijo Goldstein linesearch.\n",
+    "# In this case that linesearch causes extrapolation of the atmospheric model into invalid regions.\n",
+    "# To account for this, we explicitly assign a NewtonSolver with a different linesearch.\n",
+    "phase.nonlinear_solver = om.NewtonSolver(iprint=0, solve_subsystems=True)\n",
+    "phase.nonlinear_solver.linesearch = None\n",
+    "\n",
     "phase.add_objective('r', loc='final', scaler=-1.0)\n",
     "\n",
     "p.model.connect('size_comp.mass', 'traj.phase.parameters:m')\n",
@@ -350,7 +366,7 @@
     "#####################################################\n",
     "# Run the optimization and final explicit simulation\n",
     "#####################################################\n",
-    "dm.run_problem(p)"
+    "dm.run_problem(p, simulate=True, make_plots=True)"
    ]
   },
   {
@@ -363,60 +379,21 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "scrolled": false
+   },
    "outputs": [],
    "source": [
-    "exp_out = traj.simulate()\n",
-    "\n",
-    "rad = p.get_val('radius', units='m')[0]\n",
-    "print(f'optimal radius: {rad} m ')\n",
-    "mass = p.get_val('size_comp.mass', units='kg')[0]\n",
-    "print(f'cannonball mass: {mass} kg ')\n",
-    "angle = p.get_val('traj.phase.timeseries.gam', units='deg')[0, 0]\n",
-    "print(f'launch angle: {angle} deg')\n",
-    "max_range = p.get_val('traj.phase.timeseries.r')[-1, 0]\n",
-    "print(f'maximum range: {max_range} m')\n",
-    "\n",
-    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))\n",
-    "\n",
-    "time_imp = p.get_val('traj.phase.timeseries.time')\n",
-    "\n",
-    "time_exp = exp_out.get_val('traj.phase.timeseries.time')\n",
-    "\n",
-    "r_imp = p.get_val('traj.phase.timeseries.r')\n",
-    "\n",
-    "r_exp = exp_out.get_val('traj.phase.timeseries.r')\n",
-    "\n",
-    "h_imp = p.get_val('traj.phase.timeseries.h')\n",
-    "\n",
-    "h_exp = exp_out.get_val('traj.phase.timeseries.h')\n",
-    "\n",
-    "axes.plot(r_imp, h_imp, 'bo')\n",
-    "\n",
-    "axes.plot(r_exp, h_exp, 'r--')\n",
-    "\n",
-    "axes.set_xlabel('range (m)')\n",
-    "axes.set_ylabel('altitude (m)')\n",
-    "\n",
-    "fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))\n",
-    "states = ['r', 'h', 'v', 'gam']\n",
-    "for i, state in enumerate(states):\n",
-    "    x_imp = p.get_val(f'traj.phase.timeseries.{state}')\n",
-    "\n",
-    "    x_exp = exp_out.get_val(f'traj.phase.timeseries.{state}')\n",
-    "\n",
-    "    axes[i].set_ylabel(state)\n",
-    "\n",
-    "    axes[i].plot(time_imp, x_imp, 'bo')\n",
-    "    axes[i].plot(time_exp, x_exp, 'r--')\n",
+    "from IPython.display import IFrame\n",
     "\n",
-    "plt.show()"
+    "IFrame(src='reports/cannonball_implicit_duration/traj_results_report.html', width=700, height=1000)"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {
+    "scrolled": true,
     "tags": [
      "remove-input",
      "remove-output"
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py b/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
index f09094896..aa9110de2 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
@@ -9,6 +9,11 @@
 from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
+try:
+    from petsc4py import PETSc
+except ImportError:
+    PETSc = None
+
 from openmdao.utils.general_utils import set_pyoptsparse_opt
 
 OPT, OPTIMIZER = set_pyoptsparse_opt('SLSQP')
@@ -16,7 +21,8 @@
 
 def _make_problem(transcription='gauss-lobatto', num_segments=8, transcription_order=3,
                   compressed=True, optimizer='SLSQP', force_alloc_complex=False,
-                  solve_segments=False, y_bounds=None):
+                  solve_segments=False, y_bounds=None, phase_nonlinear_solver=None,
+                  phase_linear_solver=None):
     p = om.Problem(model=om.Group())
 
     p.driver = om.pyOptSparseDriver()
@@ -44,6 +50,12 @@ def _make_problem(transcription='gauss-lobatto', num_segments=8, transcription_o
 
     traj = dm.Trajectory()
     phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)
+
+    if phase_nonlinear_solver is not None:
+        phase.nonlinear_solver = phase_nonlinear_solver
+    if phase_linear_solver is not None:
+        phase.linear_solver = phase_linear_solver
+
     p.model.add_subsystem('traj0', traj)
     traj.add_phase('phase0', phase)
 
@@ -461,3 +473,90 @@ def test_brachistochrone_bounded_solve_segments(self):
 
             warnings.simplefilter('always')
             self.assertIn(expected_warning, [str(w.message) for w in ctx])
+
+    @unittest.skipIf(PETSc is None, 'PETSc is not available.')
+    def test_brachistochrone_solve_segments_radau_krylov_solver(self, optimizer='SLSQP'):
+        #
+        # Initialize the Problem and the optimization driver
+        #
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = optimizer
+        if optimizer == 'SNOPT':
+            p.driver.opt_settings['iSumm'] = 6
+            p.driver.opt_settings['Verify level'] = 3
+        elif optimizer == 'IPOPT':
+            p.driver.opt_settings['mu_init'] = 1e-3
+            p.driver.opt_settings['max_iter'] = 500
+            p.driver.opt_settings['print_level'] = 0
+            p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'  # for faster convergence
+            p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
+            p.driver.opt_settings['mu_strategy'] = 'monotone'
+        p.driver.declare_coloring(tol=1.0E-12)
+
+        #
+        # Create a trajectory and add a phase to it
+        #
+        traj = p.model.add_subsystem('traj0', dm.Trajectory())
+
+        # -------------------------------
+        # forward solver-based shooting
+        # -------------------------------
+        phase0 = dm.Phase(ode_class=BrachistochroneODE,
+                          transcription=dm.GaussLobatto(num_segments=10, compressed=True,
+                                                        solve_segments='forward'))
+        phase = traj.add_phase('phase0', phase0)
+
+        #
+        # Set the variables
+        #
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=False)
+
+        phase.add_state('y', fix_initial=True, fix_final=False)
+
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_control('theta', continuity=True, rate_continuity=True,
+                          units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', val=9.80665)
+
+        #
+        # Minimize time at the end of the phase
+        #
+        phase.add_objective('time', loc='final', scaler=10)
+        phase.add_boundary_constraint('x', loc='final', equals=10)
+        phase.add_boundary_constraint('y', loc='final', equals=5)
+
+        # ---------------------------------------------
+        # change linear solver for shooting
+        # ---------------------------------------------
+        # linear solver
+        phase0.linear_solver = om.PETScKrylov(assemble_jac=False, iprint=-1)
+        phase0.linear_solver.precon = om.LinearRunOnce(iprint=-1)
+
+        #
+        # Setup the Problem
+        #
+        p.setup()
+
+        #
+        # Set the initial values
+        #
+        p['traj0.phase0.t_initial'] = 0.0
+        p['traj0.phase0.t_duration'] = 2.0
+
+        p.set_val('traj0.phase0.states:x', phase.interp('x', ys=[0, 10]))
+        p.set_val('traj0.phase0.states:y', phase.interp('y', ys=[10, 5]))
+        p.set_val('traj0.phase0.states:v', phase.interp('v', ys=[0, 9.9]))
+        p.set_val('traj0.phase0.controls:theta', phase.interp('theta', ys=[5, 100.5]))
+
+        #
+        # Solve for the optimal trajectory
+        #
+        dm.run_problem(p, run_driver=True)
+
+        self.assert_results(p)
diff --git a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py
index bc607b7aa..fe61c2514 100644
--- a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py
+++ b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py
@@ -1,14 +1,10 @@
 import unittest
 
-import matplotlib.pyplot as plt
 import numpy as np
 
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 
-plt.switch_backend('Agg')
-
-
 def initial_guess(t_dur, gam0=np.pi/3):
     t = np.linspace(0, t_dur, num=100)
     g = -9.80665
@@ -206,6 +202,12 @@ def compute(self, inputs, outputs):
         # converge to the initial point.
         phase.set_duration_balance('h', val=0.0)
 
+        # In this problem, the default ArmijoGoldsteinLS has issues with extrapolating
+        # the states and causes the optimization to fail.
+        # Using the default linesearch or BoundsEnforceLS work better here.
+        phase.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, iprint=2)
+        phase.nonlinear_solver.linesearch = None
+
         phase.add_objective('r', loc='final', scaler=-1.0)
 
         p.model.connect('size_comp.mass', 'traj.phase.parameters:m')
@@ -235,60 +237,11 @@ def compute(self, inputs, outputs):
         #####################################################
         # Run the optimization and final explicit simulation
         #####################################################
-        dm.run_problem(p)
+        dm.run_problem(p, simulate=True)
 
         assert_near_equal(p.get_val('traj.phase.states:r')[-1],
                           3183.25, tolerance=1.0)
 
-        exp_out = traj.simulate()
-
-        #############################################
-        # Plot the results
-        #############################################
-        rad = p.get_val('radius', units='m')[0]
-        print(f'optimal radius: {rad} m ')
-        mass = p.get_val('size_comp.mass', units='kg')[0]
-        print(f'cannonball mass: {mass} kg ')
-        angle = p.get_val('traj.phase.timeseries.gam', units='deg')[0, 0]
-        print(f'launch angle: {angle} deg')
-        max_range = p.get_val('traj.phase.timeseries.r')[-1, 0]
-        print(f'maximum range: {max_range} m')
-
-        fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))
-
-        time_imp = p.get_val('traj.phase.timeseries.time')
-
-        time_exp = exp_out.get_val('traj.phase.timeseries.time')
-
-        r_imp = p.get_val('traj.phase.timeseries.r')
-
-        r_exp = exp_out.get_val('traj.phase.timeseries.r')
-
-        h_imp = p.get_val('traj.phase.timeseries.h')
-
-        h_exp = exp_out.get_val('traj.phase.timeseries.h')
-
-        axes.plot(r_imp, h_imp, 'bo')
-
-        axes.plot(r_exp, h_exp, 'b--')
-
-        axes.set_xlabel('range (m)')
-        axes.set_ylabel('altitude (m)')
-
-        fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))
-        states = ['r', 'h', 'v', 'gam']
-        for i, state in enumerate(states):
-            x_imp = p.get_val(f'traj.phase.timeseries.{state}')
-
-            x_exp = exp_out.get_val(f'traj.phase.timeseries.{state}')
-
-            axes[i].set_ylabel(state)
-
-            axes[i].plot(time_imp, x_imp, 'bo')
-            axes[i].plot(time_exp, x_exp, 'b--')
-
-        plt.show()
-
 
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()
diff --git a/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py b/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py
index 78c3a2efc..80fa11393 100644
--- a/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py
+++ b/dymos/examples/cannonball/doc/test_doc_two_phase_cannonball.py
@@ -1,13 +1,8 @@
 import unittest
 
-import matplotlib.pyplot as plt
-
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 
-plt.switch_backend('Agg')
-
-
 @use_tempdirs
 class TestTwoPhaseCannonballForDocs(unittest.TestCase):
 
@@ -257,75 +252,11 @@ def compute(self, inputs, outputs):
         #####################################################
         # Run the optimization and final explicit simulation
         #####################################################
-        dm.run_problem(p)
+        dm.run_problem(p, simulate=True)
 
         assert_near_equal(p.get_val('traj.descent.states:r')[-1],
                           3183.25, tolerance=1.0E-2)
 
-        # use the explicit simulation to check the final collocation solution accuracy
-        exp_out = traj.simulate()
-
-        #############################################
-        # Plot the results
-        #############################################
-        rad = p.get_val('radius', units='m')[0]
-        print(f'optimal radius: {rad} m ')
-        mass = p.get_val('size_comp.mass', units='kg')[0]
-        print(f'cannonball mass: {mass} kg ')
-        angle = p.get_val('traj.ascent.timeseries.gam', units='deg')[0, 0]
-        print(f'launch angle: {angle} deg')
-        max_range = p.get_val('traj.descent.timeseries.r')[-1, 0]
-        print(f'maximum range: {max_range} m')
-
-        fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))
-
-        time_imp = {'ascent': p.get_val('traj.ascent.timeseries.time'),
-                    'descent': p.get_val('traj.descent.timeseries.time')}
-
-        time_exp = {'ascent': exp_out.get_val('traj.ascent.timeseries.time'),
-                    'descent': exp_out.get_val('traj.descent.timeseries.time')}
-
-        r_imp = {'ascent': p.get_val('traj.ascent.timeseries.r'),
-                 'descent': p.get_val('traj.descent.timeseries.r')}
-
-        r_exp = {'ascent': exp_out.get_val('traj.ascent.timeseries.r'),
-                 'descent': exp_out.get_val('traj.descent.timeseries.r')}
-
-        h_imp = {'ascent': p.get_val('traj.ascent.timeseries.h'),
-                 'descent': p.get_val('traj.descent.timeseries.h')}
-
-        h_exp = {'ascent': exp_out.get_val('traj.ascent.timeseries.h'),
-                 'descent': exp_out.get_val('traj.descent.timeseries.h')}
-
-        axes.plot(r_imp['ascent'], h_imp['ascent'], 'bo')
-
-        axes.plot(r_imp['descent'], h_imp['descent'], 'ro')
-
-        axes.plot(r_exp['ascent'], h_exp['ascent'], 'b--')
-
-        axes.plot(r_exp['descent'], h_exp['descent'], 'r--')
-
-        axes.set_xlabel('range (m)')
-        axes.set_ylabel('altitude (m)')
-
-        fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))
-        states = ['r', 'h', 'v', 'gam']
-        for i, state in enumerate(states):
-            x_imp = {'ascent': p.get_val(f'traj.ascent.timeseries.{state}'),
-                     'descent': p.get_val(f'traj.descent.timeseries.{state}')}
-
-            x_exp = {'ascent': exp_out.get_val(f'traj.ascent.timeseries.{state}'),
-                     'descent': exp_out.get_val(f'traj.descent.timeseries.{state}')}
-
-            axes[i].set_ylabel(state)
-
-            axes[i].plot(time_imp['ascent'], x_imp['ascent'], 'bo')
-            axes[i].plot(time_imp['descent'], x_imp['descent'], 'ro')
-            axes[i].plot(time_exp['ascent'], x_exp['ascent'], 'b--')
-            axes[i].plot(time_exp['descent'], x_exp['descent'], 'r--')
-
-        plt.show()
-
 
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()
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 c6026d709..75621f3b1 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
@@ -10,7 +10,7 @@
 def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=False,
               connected=False, default_nonlinear_solver=None, default_linear_solver=None):
     """
-    Build a traejctory for the finite burn orbit raise problem.
+    Build a trajectory for the finite burn orbit raise problem.
 
     Parameters
     ----------
@@ -31,8 +31,13 @@ def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=F
     t = {'gauss-lobatto': dm.GaussLobatto(num_segments=5, order=transcription_order, compressed=compressed),
          'radau': dm.Radau(num_segments=5, order=transcription_order, compressed=compressed)}
 
-    traj = dm.Trajectory(default_nonlinear_solver=default_nonlinear_solver,
-                         default_linear_solver=default_linear_solver)
+    traj = dm.Trajectory()
+
+    if default_nonlinear_solver is not None:
+        traj.phases.nonlinear_solver = default_nonlinear_solver
+
+    if default_linear_solver is not None:
+        traj.phases.linear_solver = default_linear_solver
 
     traj.add_parameter('c', opt=False, val=1.5, units='DU/TU',
                        targets={'burn1': ['c'], 'burn2': ['c']})
diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py b/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py
index c3d11c8d2..ad55821bc 100644
--- a/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py
+++ b/dymos/examples/finite_burn_orbit_raise/test/test_multi_phase_restart.py
@@ -114,7 +114,6 @@ def test_ex_two_burn_orbit_raise_connected(self):
                                          show_output=False, restart='dymos_solution.db',
                                          connected=True, solution_record_file='dymos_solution2.db',
                                          simulation_record_file='dymos_simulation2.db')
-#
 
         case2 = om.CaseReader('dymos_solution2.db').get_case('final')
         sim_case2 = om.CaseReader('dymos_simulation2.db').get_case('final')
@@ -127,41 +126,30 @@ def test_ex_two_burn_orbit_raise_connected(self):
     def test_restart_from_solution_radau_to_connected(self):
         optimizer = 'IPOPT'
 
-        expected_warnings = \
-            [(om.OpenMDAOWarning,
-              "'traj' <class Trajectory>: Setting phases.nonlinear_solver to `om.NonlinearBlockJac(iprint=0)`.\n"
-              "Connected phases in parallel require a non-default nonlinear solver.\n"
-              "Use traj.options[\'default_nonlinear_solver\'] to explicitly set the solver."),
-             (om.OpenMDAOWarning,
-              "'traj' <class Trajectory>: Setting phases.linear_solver to `om.PETScKrylov()`.\n"
-              "Connected phases in parallel require a non-default linear solver.\n"
-              "Use traj.options[\'default_linear_solver\'] to explicitly set the solver.")]
-
-        with assert_warnings(expected_warnings):
-            p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3,
-                                             compressed=False, optimizer=optimizer,
-                                             show_output=False, connected=True)
+        p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3,
+                                         compressed=False, optimizer=optimizer,
+                                         show_output=False, connected=True)
 
-            if p.model.traj.phases.burn2 in p.model.traj.phases._subsystems_myproc:
-                assert_near_equal(p.get_val('traj.burn2.states:deltav')[0], 0.3995,
-                                  tolerance=4.0E-3)
+        if p.model.traj.phases.burn2 in p.model.traj.phases._subsystems_myproc:
+            assert_near_equal(p.get_val('traj.burn2.states:deltav')[0], 0.3995,
+                              tolerance=4.0E-3)
 
-            case1 = om.CaseReader('dymos_solution.db').get_case('final')
-            sim_case1 = om.CaseReader('dymos_simulation.db').get_case('final')
+        case1 = om.CaseReader('dymos_solution.db').get_case('final')
+        sim_case1 = om.CaseReader('dymos_simulation.db').get_case('final')
 
-            # Run again without an actual optimizer
-            p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3,
-                                             compressed=False, optimizer=optimizer, run_driver=False,
-                                             show_output=False, restart='dymos_solution.db',
-                                             connected=True, solution_record_file='dymos_solution2.db',
-                                             simulation_record_file='dymos_simulation2.db')
+        # Run again without an actual optimizer
+        p = two_burn_orbit_raise_problem(transcription='radau', transcription_order=3,
+                                         compressed=False, optimizer=optimizer, run_driver=False,
+                                         show_output=False, restart='dymos_solution.db',
+                                         connected=True, solution_record_file='dymos_solution2.db',
+                                         simulation_record_file='dymos_simulation2.db')
 
-            case2 = om.CaseReader('dymos_solution2.db').get_case('final')
-            sim_case2 = om.CaseReader('dymos_simulation2.db').get_case('final')
+        case2 = om.CaseReader('dymos_solution2.db').get_case('final')
+        sim_case2 = om.CaseReader('dymos_simulation2.db').get_case('final')
 
-            # Verify that the second case has the same inputs and outputs
-            assert_cases_equal(case1, case2, tol=1.0E-8)
-            assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8)
+        # Verify that the second case has the same inputs and outputs
+        assert_cases_equal(case1, case2, tol=1.0E-8)
+        assert_cases_equal(sim_case1, sim_case2, tol=1.0E-8)
 
 
 if __name__ == '__main__':  # pragma: no cover
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 73a5f8b16..d036fafcb 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -114,6 +114,8 @@ def initialize(self):
                              desc='Keyword arguments provided when initializing the ODE System')
         self.options.declare('transcription', types=TranscriptionBase,
                              desc='Transcription technique of the optimal control problem.')
+        self.options.declare('auto_solvers', types=bool, default=True,
+                             desc='If True, attempt to automatically assign solvers if necessary.')
 
     def add_state(self, name, units=_unspecified, shape=_unspecified,
                   rate_source=_unspecified, targets=_unspecified,
diff --git a/dymos/trajectory/test/test_trajectory.py b/dymos/trajectory/test/test_trajectory.py
index 45ef246d1..c65d9080a 100644
--- a/dymos/trajectory/test/test_trajectory.py
+++ b/dymos/trajectory/test/test_trajectory.py
@@ -1779,7 +1779,7 @@ def compute(self, inputs, outputs):
         p.model.add_design_var('radius', lower=0.01, upper=0.10,
                                ref0=0.01, ref=0.10, units='m')
 
-        traj = p.model.add_subsystem('traj', dm.Trajectory())
+        traj = p.model.add_subsystem('traj', dm.Trajectory(parallel_phases=True))
 
         transcription = dm.Radau(num_segments=5, order=3, compressed=True)
 
@@ -1856,9 +1856,6 @@ def compute(self, inputs, outputs):
         p.model.connect('size_comp.mass', 'traj.parameters:m')
         p.model.connect('size_comp.S', 'traj.parameters:S')
 
-        # A linear solver at the top level can improve performance.
-        p.model.linear_solver = om.DirectSolver()
-
         # Finish Problem Setup
         p.setup()
 
diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py
index 31cce9b7d..55886fe79 100644
--- a/dymos/trajectory/trajectory.py
+++ b/dymos/trajectory/trajectory.py
@@ -61,7 +61,7 @@ def __init__(self, **kwargs):
         self._has_connected_phases = False
 
         self.parameter_options = {}
-        self.phases = om.ParallelGroup()
+        self.phases = om.ParallelGroup() if self.options['parallel_phases'] else om.Group()
 
     def initialize(self):
         """
@@ -69,17 +69,12 @@ def initialize(self):
         """
         self.options.declare('sim_mode', types=bool, default=False,
                              desc='Used internally by Dymos when invoking simulate on a trajectory')
-        self.options.declare('default_nonlinear_solver',
-                             types=(om.NonlinearBlockJac, om.NewtonSolver, om.BroydenSolver),
-                             default=None, allow_none=True, recordable=False,
-                             desc='A nonlinear solver to be used when Phases are connected but being '
-                                  'run in parallel. If not specified, Dymos will automatically use '
-                                  'NonlinearBlockJac in this situation.')
-        self.options.declare('default_linear_solver', default=None, allow_none=True, recordable=False,
-                             types=(om.PETScKrylov, om.LinearBlockJac, om.NonlinearBlockGS),
-                             desc='A linear solver to be used when Phases are connected but being '
-                                  'run in parallel. If not specified, Dymos will automatically use '
-                                  'PETScKrylov in this situation.')
+        self.options.declare('parallel_phases',
+                             default=True, types=bool,
+                             desc='If True, the top-level container of all phases will be a ParallelGroup,'
+                                  'otherwise it will be a standard OpenMDAO Group.')
+        self.options.declare('auto_solvers', types=bool, default=True,
+                             desc='If True, attempt to automatically assign solvers if necessary.')
 
     def add_phase(self, name, phase, **kwargs):
         """
@@ -1018,27 +1013,30 @@ def _configure_solvers(self):
         These solvers can be changed through the
         'default_nonlinear_solver' and 'default_linear_solver' options.
         """
-        if self._has_connected_phases and MPI and self.comm.size > 1:
+        if not self.options['auto_solvers']:
+            return
 
+        warn = False
+        if self._has_connected_phases and isinstance(self.phases, om.ParallelGroup) and self.comm.size > 1:
+            msg = (f'{self.msginfo}\n'
+                   f'  Non-default solvers are required\n'
+                   f'    connected phases in parallel: True\n')
             if isinstance(self.phases.nonlinear_solver, om.NonlinearRunOnce):
-                if self.options['default_nonlinear_solver'] is None:
-                    msg = f'{self.msginfo}: Setting phases.nonlinear_solver to `om.NonlinearBlockJac(iprint=0)`.\n' \
-                          f'Connected phases in parallel require a non-default nonlinear solver.\n' \
-                          f'Use {self.pathname}.options[\'default_nonlinear_solver\'] to explicitly set the solver.'
-                    om.issue_warning(msg)
-                    self.phases.nonlinear_solver = om.NonlinearBlockJac(iprint=0)
-                else:
-                    self.phases.nonlinear_solver = self.options['default_nonlinear_solver']
+                warn = True
+                msg += (f'  Setting \n'
+                        f'    {self.pathname}.phases.nonlinear_solver = om.NonlinearBlockJac(iprint=0)\n'
+                        f'  Explicitly set {self.pathname}.phases.nonlinear_solver to override.\n')
+                self.phases.nonlinear_solver = om.NonlinearBlockJac(iprint=0)
 
             if isinstance(self.phases.linear_solver, om.LinearRunOnce):
-                if self.options['default_linear_solver'] is None:
-                    msg = f'{self.msginfo}: Setting phases.linear_solver to `om.PETScKrylov()`.\n' \
-                          f'Connected phases in parallel require a non-default linear solver.\n' \
-                          f'Use {self.pathname}.options[\'default_linear_solver\'] to explicitly set the solver.'
-                    om.issue_warning(msg)
-                    self.phases.linear_solver = om.PETScKrylov()
-                else:
-                    self.phases.linear_solver = self.options['default_linear_solver']
+                warn = True
+                msg += (f'  Setting\n'
+                        f'    {self.pathname}.phases.linear_solver = om.PETScKrylov(iprint=0)\n'
+                        f'  Explicitly set {self.pathname}.phases.linear_solver to override.\n')
+                self.phases.linear_solver = om.PETScKrylov(iprint=0)
+
+            if warn:
+                om.issue_warning(msg)
 
     def configure(self):
         """
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index 1c7d28246..e51342eb8 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -10,7 +10,6 @@
 from ..transcription_base import TranscriptionBase
 from ..grid_data import BirkhoffGrid, GaussLobattoGrid, RadauGrid, UniformGrid
 from .ode_integration_comp import ODEIntegrationComp
-from ..._options import options as dymos_options
 from ...utils.misc import get_rate_units, CoerceDesvar
 from ...utils.indexing import get_src_indices_by_row
 from ...utils.introspection import get_promoted_vars, get_source_metadata, get_targets, _get_targets_metadata
@@ -675,7 +674,7 @@ def setup_solvers(self, phase):
         """
         pass
 
-    def configure_solvers(self, phase):
+    def configure_solvers(self, phase, requires_solvers=None):
         """
         Setup the solvers for this transcription.
 
@@ -683,8 +682,10 @@ def configure_solvers(self, phase):
         ----------
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
+        requires_solvers : None
+            Required to extend TranscriptionBase.configure_solvers but not used.
         """
-        pass
+        super().configure_solvers(phase, requires_solvers=None)
 
     def get_parameter_connections(self, name, phase):
         """
diff --git a/dymos/transcriptions/pseudospectral/birkhoff.py b/dymos/transcriptions/pseudospectral/birkhoff.py
index df914ead0..9840ab150 100644
--- a/dymos/transcriptions/pseudospectral/birkhoff.py
+++ b/dymos/transcriptions/pseudospectral/birkhoff.py
@@ -343,7 +343,7 @@ def setup_solvers(self, phase):
         """
         pass
 
-    def configure_solvers(self, phase):
+    def configure_solvers(self, phase, requires_solvers=None):
         """
         Configure the solvers.
 
@@ -351,35 +351,17 @@ def configure_solvers(self, phase):
         ----------
         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 ode_iter_group._implicit_outputs or self._implicit_duration:
-            # Only override the solvers if the user hasn't set them to something else.
-            if isinstance(phase.nonlinear_solver, om.NonlinearRunOnce):
-                msg = f'{phase.msginfo}: Phase requires non-default nonlinear solver due to the use of\n' \
-                      f'solve_segments or the use of set_duration_balance.\n' \
-                      f'Setting nonlinear solver to ' \
-                      f'om.NewtonSolver(solve_subsystems=True, maxiter=100, iprint=2, stall_limit=3)\n' \
-                      f'by default.' \
-                      f'To avoid this warning explicitly assign a nonlinear solver to this phase.'
-
-                om.issue_warning(msg)
-                newton = phase.nonlinear_solver = om.NewtonSolver()
-                newton.options['solve_subsystems'] = True
-                newton.options['maxiter'] = 100
-                newton.options['iprint'] = 2
-                newton.options['stall_limit'] = 3
-                newton.linesearch = om.BoundsEnforceLS()
-
-            if isinstance(phase.linear_solver, om.LinearRunOnce):
-                msg = f'{phase.msginfo}: Phase requires non-default linear solver due to the use of\n' \
-                      f'solve_segments or the use of set_duration_balance.\n' \
-                      f'Setting linear solver to om.DirectSolver() by default.' \
-                      f'To avoid this warning explicitly assign a linear solver to this phase.'
-
-                om.issue_warning(msg)
-                phase.linear_solver = om.DirectSolver()
+        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):
         """
diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
index c9b9f9625..411440065 100644
--- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py
+++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
@@ -20,6 +20,12 @@ class PseudospectralBase(TranscriptionBase):
     **kwargs : dict
         Dictionary of optional arguments.
     """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+        self.any_solved_segs = False
+        self.any_connected_opt_segs = False
+
     def initialize(self):
         """
         Declare transcription options.
@@ -80,8 +86,6 @@ def setup_states(self, phase):
         """
         grid_data = self.grid_data
 
-        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:
@@ -95,7 +99,6 @@ def setup_states(self, phase):
         if self.any_solved_segs or self.any_connected_opt_segs:
             indep = StateIndependentsComp(grid_data=grid_data,
                                           state_options=phase.state_options)
-            indep.linear_solver = om.DirectSolver()
         else:
             indep = om.IndepVarComp()
 
@@ -495,7 +498,7 @@ def setup_solvers(self, phase):
         """
         pass
 
-    def configure_solvers(self, phase):
+    def configure_solvers(self, phase, requires_solvers=None):
         """
         Configure the solvers.
 
@@ -503,19 +506,15 @@ def configure_solvers(self, phase):
         ----------
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
-        """
-        if self.any_solved_segs or self._implicit_duration:
-            # Only override the solvers if the user hasn't set them to something else.
-            if isinstance(phase.nonlinear_solver, om.NonlinearRunOnce):
-                newton = phase.nonlinear_solver = om.NewtonSolver()
-                newton.options['solve_subsystems'] = True
-                newton.options['maxiter'] = 100
-                newton.options['iprint'] = 2
-                newton.options['stall_limit'] = 3
-                newton.linesearch = om.BoundsEnforceLS()
-
-            if isinstance(phase.linear_solver, om.LinearRunOnce):
-                phase.linear_solver = om.DirectSolver()
+        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.
+        """
+        req_solvers = {'solved segments': self.any_solved_segs,
+                       'input initial': self.any_connected_opt_segs}
+        if requires_solvers is not None:
+            req_solvers.update(requires_solvers)
+        super().configure_solvers(phase, requires_solvers=req_solvers)
 
     def setup_timeseries_outputs(self, phase):
         """
diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index db64a0a05..077fb1901 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -399,17 +399,53 @@ def setup_solvers(self, phase):
         """
         raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement method setup_solvers.')
 
-    def configure_solvers(self, phase):
+    def configure_solvers(self, phase, requires_solvers=None):
         """
-        Configure the solvers for this transcription.
+        Configure the solvers.
 
         Parameters
         ----------
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
-        """
-        raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement method '
-                                  f'configure_solvers.')
+        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.
+        """
+        if not phase.options['auto_solvers']:
+            return
+
+        req_solvers = {'implicit duration': self._implicit_duration}
+
+        if requires_solvers is not None:
+            req_solvers.update(requires_solvers)
+
+        reasons_txt = '\n'.join(f'    {key}: {val}' for key, val in req_solvers.items())
+        warn = False
+
+        if any(req_solvers.values()):
+            msg = (f'{phase.msginfo}\n'
+                   f'  Non-default solvers are required\n'
+                   f'{reasons_txt}\n')
+            if isinstance(phase.nonlinear_solver, om.NonlinearRunOnce):
+                msg += (f'  Setting `{phase.pathname}.nonlinear_solver = om.NewtonSolver(iprint=0, '
+                        f'solve_subsystems=True, maxiter=1000, stall_limit=3)`\n'
+                        f'  Explicitly set {phase.pathname}.nonlinear_solver to override.\n')
+                warn = True
+                phase.nonlinear_solver = om.NewtonSolver(iprint=0)
+                phase.nonlinear_solver.options['solve_subsystems'] = True
+                phase.nonlinear_solver.options['maxiter'] = 1000
+                phase.nonlinear_solver.options['stall_limit'] = 3
+                phase.nonlinear_solver.linesearch = om.ArmijoGoldsteinLS()
+
+            if isinstance(phase.linear_solver, om.LinearRunOnce):
+                warn = True
+                msg += (f'  Setting `{phase.pathname}.linear_solver = om.DirectSolver(iprint=2)`\n'
+                        f'  Explicitly set {phase.pathname}.linear_solver to override.\n')
+                phase.linear_solver = om.DirectSolver(iprint=0)
+
+            if warn:
+                msg += f'  Set `{phase.pathname}.options["auto_solvers"] = False` to disable this behavior.'
+                om.issue_warning(msg)
 
     def setup_timeseries_outputs(self, phase):
         """