diff --git a/.gitignore b/.gitignore
index 650d350f3..c2a245e63 100644
--- a/.gitignore
+++ b/.gitignore
@@ -81,10 +81,15 @@ fort.6
 # mkdocs
 /mkdocs/site
 
-# Examples
-dymos/examples/*/*/_output/
-docs/examples/*/*.db
-docs/dymos_book/examples/*/coloring_files/
+# Docs
+docs/dymos_book/_build
+docs/dymos_book/_srcdocs
+
+# Dymos and OpenMDAO files
+**/*.db
+**/*.sql
+**/coloring_files
+**/reports
 
 # Misc files
 docs/faq/*/*.pkl
diff --git a/docs/dymos_book/api/phase_api.ipynb b/docs/dymos_book/api/phase_api.ipynb
index 5eeaf233f..01ca8d272 100644
--- a/docs/dymos_book/api/phase_api.ipynb
+++ b/docs/dymos_book/api/phase_api.ipynb
@@ -248,6 +248,12 @@
     "        :noindex:\n",
     "```\n",
     "\n",
+    "## duplicate\n",
+    "```{eval-rst}\n",
+    "    .. automethod:: dymos.Phase.duplicate\n",
+    "        :noindex:\n",
+    "```\n",
+    "\n",
     "## set_refine_options\n",
     "```{eval-rst}\n",
     "    .. automethod:: dymos.Phase.set_refine_options\n",
diff --git a/docs/dymos_book/examples/multi_phase_cannonball/multi_phase_cannonball.ipynb b/docs/dymos_book/examples/multi_phase_cannonball/multi_phase_cannonball.ipynb
index 10c987b17..b746b2b2d 100644
--- a/docs/dymos_book/examples/multi_phase_cannonball/multi_phase_cannonball.ipynb
+++ b/docs/dymos_book/examples/multi_phase_cannonball/multi_phase_cannonball.ipynb
@@ -350,24 +350,16 @@
     "\n",
     "# Second Phase (descent)\n",
     "transcription = dm.GaussLobatto(num_segments=5, order=3, compressed=True)\n",
-    "descent = dm.Phase(ode_class=CannonballODE, transcription=transcription)\n",
+    "# descent = dm.Phase(ode_class=CannonballODE, transcription=transcription)\n",
+    "descent = ascent.duplicate(transcription=transcription)\n",
     "\n",
     "traj.add_phase('descent', descent)\n",
     "\n",
-    "# All initial states and time are free, since\n",
-    "#    they will be linked to the final states of ascent.\n",
-    "# Final altitude is fixed, because we will set\n",
-    "#    it to zero so that the phase ends at ground impact)\n",
-    "descent.set_time_options(initial_bounds=(.5, 100), duration_bounds=(.5, 100),\n",
-    "                         duration_ref=100, units='s')\n",
-    "descent.add_state('r')\n",
-    "descent.add_state('h', fix_initial=False, fix_final=True)\n",
-    "descent.add_state('gam', fix_initial=False, fix_final=False)\n",
-    "descent.add_state('v', fix_initial=False, fix_final=False)\n",
-    "\n",
-    "descent.add_parameter('S', units='m**2', static_target=True)\n",
-    "descent.add_parameter('m', units='kg', static_target=True)\n",
-    "\n",
+    "# Because we copied the descent phase\n",
+    "# - The 'fix_initial' option for time was set to False\n",
+    "# - All state 'fix_initial' and 'fix_final' options are set to False.\n",
+    "# - We only need to fix the final value of altitude so the descent phase ends at ground impact.\n",
+    "descent.set_state_options('h', fix_final=True)\n",
     "descent.add_objective('r', loc='final', scaler=-1.0)\n",
     "\n",
     "# Add internally-managed design parameters to the trajectory.\n",
@@ -554,7 +546,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.7"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,
diff --git a/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py b/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py
index 17ecf9157..db453d5c7 100644
--- a/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py
+++ b/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py
@@ -96,7 +96,7 @@ def setup(self):
         rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]
         self.rho_interp = interp1d(np.array(alt_data, dtype=complex),
                                    np.array(rho_data, dtype=complex),
-                                   kind='linear')
+                                   kind='linear', bounds_error=False, fill_value='extrapolate')
 
     def compute(self, inputs, outputs):
 
diff --git a/dymos/phase/analytic_phase.py b/dymos/phase/analytic_phase.py
index 3d6212c93..eb3b416c6 100644
--- a/dymos/phase/analytic_phase.py
+++ b/dymos/phase/analytic_phase.py
@@ -1,3 +1,5 @@
+from copy import deepcopy
+
 import openmdao.api as om
 
 from .phase import Phase
@@ -472,3 +474,67 @@ def set_simulate_options(self, *args, **kwargs):
             Simulation cannot be performed on AnalyticPhase.
         """
         raise NotImplementedError('Method set_simulate_options is not available for AnalyticPhase.')
+
+    def duplicate(self, num_nodes=None, boundary_constraints=False, path_constraints=False, objectives=False,
+                  fix_initial_time=False):
+        """
+        Create a copy of this phase where most options and attributes are deep copies of those in the original.
+
+        By default, a deepcopy of the transcription in the original phase is used.
+        Boundary constraints, path constraints, and objectives are _NOT_ copied by default, but the user may opt to do so.
+        By default, initial time is not fixed, nor are the initial or final state values.
+        These also can be overridden with the appropriate arguments.
+
+        Parameters
+        ----------
+        num_nodes : int or None
+            The number of nodes to use in the new phase, or None if it should use the same
+            number as the phase being duplicated.
+        boundary_constraints : bool
+            If True, retain all boundary constraints from the phase to be copied.
+        path_constraints : bool
+            If True, retain all path constraints from the phase to be copied.
+        objectives : bool
+            If True, retain all objectives from the phase to be copied.
+        fix_initial_time : bool
+            If True, fix the initial time of the returned phase.
+
+        Returns
+        -------
+        AnalyticPhase
+            The new phase created by duplicating this one.
+        """
+        nn = num_nodes if num_nodes is not None else self.options['num_nodes']
+        ode_class = self.options['ode_class']
+        ode_init_kwargs = self.options['ode_init_kwargs']
+        auto_solvers = self.options['auto_solvers']
+
+        p = AnalyticPhase(num_nodes=nn, ode_class=ode_class, ode_init_kwargs=ode_init_kwargs,
+                          auto_solvers=auto_solvers)
+
+        p.time_options.update(deepcopy(self.time_options))
+        p.time_options['fix_initial'] = fix_initial_time
+
+        for state_name, state_options in self.state_options.items():
+            p.state_options[state_name] = deepcopy(state_options)
+
+        for param_name, param_options in self.parameter_options.items():
+            p.parameter_options[param_name] = deepcopy(param_options)
+
+        p._timeseries = deepcopy(self._timeseries)
+
+        p.refine_options = deepcopy(self.refine_options)
+        p.simulate_options = deepcopy(self.simulate_options)
+        p.timeseries_options = deepcopy(self.timeseries_options)
+
+        if boundary_constraints:
+            p._initial_boundary_constraints = deepcopy(self._initial_boundary_constraints)
+            p._final_boundary_constraints = deepcopy(self._final_boundary_constraints)
+
+        if path_constraints:
+            p._path_constraints = deepcopy(self._path_constraints)
+
+        if objectives:
+            p._objectives = deepcopy(self._objectives)
+
+        return p
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index d036fafcb..ad3ff612d 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1,5 +1,7 @@
 from collections.abc import Iterable, Callable, Sequence
+from copy import deepcopy
 import inspect
+from os import path
 import warnings
 
 import numpy as np
@@ -21,7 +23,7 @@
     TimeseriesOutputOptionsDictionary, PhaseTimeseriesOptionsDictionary
 
 from ..transcriptions.transcription_base import TranscriptionBase
-from ..transcriptions.grid_data import GaussLobattoGrid, RadauGrid, UniformGrid
+from ..transcriptions.grid_data import GaussLobattoGrid, RadauGrid, UniformGrid, BirkhoffGrid
 from ..transcriptions import ExplicitShooting, GaussLobatto, Radau
 from ..utils.indexing import get_constraint_flat_idxs
 from ..utils.introspection import configure_time_introspection, _configure_constraint_introspection, \
@@ -60,13 +62,18 @@ class Phase(om.Group):
     def __init__(self, from_phase=None, **kwargs):
         _kwargs = kwargs.copy()
 
+        if from_phase is not None:
+            raise RuntimeError('Instantiating a phase from another using the `from_phase`'
+                               ' argument is no longer avaiable.'
+                               ' Use the <phase>.duplicate() or <phase>.get_simulation_phase() methods.')
+
         # These are the options which will be set at setup time.
         # Prior to setup, the options are placed into the user_*_options dictionaries.
-        self.time_options = TimeOptionsDictionary()
-        self.state_options = {}
-        self.control_options = {}
-        self.polynomial_control_options = {}
-        self.parameter_options = {}
+        # self.time_options = TimeOptionsDictionary()
+        # self.state_options = {}
+        # self.control_options = {}
+        # self.polynomial_control_options = {}
+        # self.parameter_options = {}
         self.refine_options = GridRefinementOptionsDictionary()
         self.simulate_options = SimulateOptionsDictionary()
         self.timeseries_ec_vars = {}
@@ -74,34 +81,104 @@ def __init__(self, from_phase=None, **kwargs):
 
         # Dictionaries of variable options that are set by the user via the API
         # These will be applied over any defaults specified by decorators on the ODE
-        if from_phase is None:
-            self._initial_boundary_constraints = []
-            self._final_boundary_constraints = []
-            self._path_constraints = []
-            self._timeseries = {'timeseries': {'transcription': None,
-                                               'subset': 'all',
-                                               'outputs': {}}}
-            self._objectives = {}
-        else:
-            self.time_options.update(from_phase.time_options)
-            self.state_options = from_phase.state_options.copy()
-            self.control_options = from_phase.control_options.copy()
-            self.polynomial_control_options = from_phase.polynomial_control_options.copy()
-            self.parameter_options = from_phase.parameter_options.copy()
+        self._initial_boundary_constraints = []
+        self._final_boundary_constraints = []
+        self._path_constraints = []
+        self._timeseries = {'timeseries': {'transcription': None,
+                                           'subset': 'all',
+                                           'outputs': {}}}
+        self._objectives = {}
+
+        super(Phase, self).__init__(**_kwargs)
 
-            self.refine_options.update(from_phase.refine_options)
-            self.simulate_options.update(from_phase.simulate_options)
+    def duplicate(self, transcription=None, boundary_constraints=False, path_constraints=False, objectives=False,
+                  fix_initial_time=False, fix_initial_states=None, fix_final_states=None):
+        """
+        Create a copy of this phase where most options and attributes are deep copies of those in the original.
 
-            self._initial_boundary_constraints = from_phase._initial_boundary_constraints.copy()
-            self._final_boundary_constraints = from_phase._final_boundary_constraints.copy()
-            self._path_constraints = from_phase._path_constraints.copy()
-            self._timeseries = from_phase._timeseries.copy()
-            self._objectives = from_phase._objectives.copy()
+        By default, a deepcopy of the transcription in the original phase is used.
+        Boundary constraints, path constraints, and objectives are _NOT_ copied by default, but the user may opt to do so.
+        By default, initial time is not fixed, nor are the initial or final state values.
+        These also can be overridden with the appropriate arguments.
 
-            _kwargs['ode_class'] = from_phase.options['ode_class']
-            _kwargs['ode_init_kwargs'] = from_phase.options['ode_init_kwargs']
+        Parameters
+        ----------
+        transcription : TranscriptionBase or None
+            If given, use the specified transcription for the new phase.
+            If None, use a copy of the transcription of the phase to be copied.
+        boundary_constraints : bool
+            If True, retain all boundary constraints from the phase to be copied.
+        path_constraints : bool
+            If True, retain all path constraints from the phase to be copied.
+        objectives : bool
+            If True, retain all objectives from the phase to be copied.
+        fix_initial_time : bool
+            If True, fix the initial time of the returned phase.
+        fix_initial_states : Sequence of str or None
+            If given, set fix_initial=True for the given state names. Otherwise, all states will have fix_initial=False.
+        fix_final_states : Sequence of str or None
+            If given, set fix_final=True for the given state names. Otherwise, all states will have fix_final=False.
 
-        super(Phase, self).__init__(**_kwargs)
+        Returns
+        -------
+        Phase
+            The new phase created by duplicating this one.
+        """
+        t = deepcopy(self.options['transcription']) if transcription is None else transcription
+        ode_class = self.options['ode_class']
+        ode_init_kwargs = self.options['ode_init_kwargs']
+        auto_solvers = self.options['auto_solvers']
+
+        p = Phase(transcription=t, ode_class=ode_class, ode_init_kwargs=ode_init_kwargs,
+                  auto_solvers=auto_solvers)
+
+        # First copy over the variable information as is
+        p.time_options.update(deepcopy(self.time_options))
+
+        for state_name, state_options in self.state_options.items():
+            p.state_options[state_name] = deepcopy(state_options)
+
+        for param_name, param_options in self.parameter_options.items():
+            p.parameter_options[param_name] = deepcopy(param_options)
+
+        for control_name, control_options in self.control_options.items():
+            p.control_options[control_name] = deepcopy(control_options)
+
+        for pc_name, pc_options in self.polynomial_control_options.items():
+            p.polynomial_control_options[pc_name] = deepcopy(pc_options)
+
+        p.time_options['fix_initial'] = fix_initial_time
+
+        _fis = [] if fix_initial_states is None else fix_initial_states
+        _ffs = [] if fix_final_states is None else fix_final_states
+
+        for state_name, state_options in p.state_options.items():
+            if state_name in _fis:
+                state_options['fix_initial'] = True
+            else:
+                state_options['fix_initial'] = False
+
+            if state_name in _ffs:
+                state_options['fix_final'] = True
+            else:
+                state_options['fix_final'] = False
+
+        p._timeseries = deepcopy(self._timeseries)
+        p.refine_options = deepcopy(self.refine_options)
+        p.simulate_options = deepcopy(self.simulate_options)
+        p.timeseries_options = deepcopy(self.timeseries_options)
+
+        if boundary_constraints:
+            p._initial_boundary_constraints = deepcopy(self._initial_boundary_constraints)
+            p._final_boundary_constraints = deepcopy(self._final_boundary_constraints)
+
+        if path_constraints:
+            p._path_constraints = deepcopy(self._path_constraints)
+
+        if objectives:
+            p._objectives = deepcopy(self._objectives)
+
+        return p
 
     def initialize(self):
         """
@@ -116,6 +193,40 @@ def initialize(self):
                              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.')
+        self.options.declare('time_options', types=TimeOptionsDictionary, default=TimeOptionsDictionary(),
+                             desc='Options for time in this phase.')
+        self.options.declare('state_options', types=dict, default={},
+                             desc='Options for each state in this phase.')
+        self.options.declare('parameter_options', types=dict, default={},
+                             desc='Options for each parameter in this phase.')
+        self.options.declare('control_options', types=dict, default={},
+                             desc='Options for each control in this phase.')
+        self.options.declare('polynomial_control_options', types=dict, default={},
+                             desc='Options for each polynomial control in this phase.')
+
+    @property
+    def time_options(self):
+        return self.options['time_options']
+
+    @time_options.setter
+    def time_options(self, tod: TimeOptionsDictionary):
+        self.options['time_options'] = tod
+
+    @property
+    def state_options(self):
+        return self.options['state_options']
+
+    @property
+    def parameter_options(self):
+        return self.options['parameter_options']
+
+    @property
+    def control_options(self):
+        return self.options['control_options']
+
+    @property
+    def polynomial_control_options(self):
+        return self.options['polynomial_control_options']
 
     def add_state(self, name, units=_unspecified, shape=_unspecified,
                   rate_source=_unspecified, targets=_unspecified,
@@ -2254,9 +2365,92 @@ def get_simulation_phase(self, times_per_seg=_unspecified, method=_unspecified,
         """
         from .simulation_phase import SimulationPhase
 
-        sim_phase = SimulationPhase(from_phase=self, times_per_seg=times_per_seg, method=method,
-                                    atol=atol, rtol=rtol, first_step=first_step, max_step=max_step,
-                                    reports=reports)
+        # t = deepcopy(self.options['transcription']) if transcription is None else transcription
+        ode_class = self.options['ode_class']
+        ode_init_kwargs = self.options['ode_init_kwargs']
+        auto_solvers = self.options['auto_solvers']
+
+        self_tx = self.options['transcription']
+        num_seg = self_tx.grid_data.num_segments
+        seg_order = self_tx.grid_data.transcription_order
+        seg_ends = self_tx.grid_data.segment_ends
+        compressed = self_tx.grid_data.compressed
+
+        sim_options = self.simulate_options
+
+        _method = method if method is not _unspecified else sim_options['method']
+        _atol = atol if atol is not _unspecified else sim_options['atol']
+        _rtol = rtol if rtol is not _unspecified else sim_options['rtol']
+        _first_step = first_step if first_step is not _unspecified else sim_options['first_step']
+        _max_step = max_step if max_step is not _unspecified else sim_options['max_step']
+        _times_per_seg = times_per_seg if times_per_seg is not _unspecified else sim_options['times_per_seg']
+
+        if isinstance(self_tx, GaussLobatto):
+            grid = GaussLobattoGrid(num_segments=num_seg, nodes_per_seg=seg_order, segment_ends=seg_ends,
+                                    compressed=compressed)
+        elif isinstance(self_tx, Radau):
+            grid = RadauGrid(num_segments=num_seg, nodes_per_seg=seg_order + 1, segment_ends=seg_ends,
+                             compressed=compressed)
+        elif isinstance(self_tx.grid_data, GaussLobattoGrid) or \
+                isinstance(self_tx.grid_data, RadauGrid) or isinstance(self_tx.grid_data, BirkhoffGrid):
+            grid = self_tx.grid_data
+        else:
+            raise RuntimeError(f'Unexpected grid class for {phase_tx.grid_data}. Only phases with GaussLobatto '
+                               f'or Radau grids can be simulated.')
+
+        if _times_per_seg is None:
+            output_grid = None
+        else:
+            output_grid = UniformGrid(num_segments=num_seg, nodes_per_seg=_times_per_seg, segment_ends=seg_ends,
+                                      compressed=compressed)
+
+        tx = ExplicitShooting(propagate_derivs=False,
+                              subprob_reports=reports,
+                              grid=grid,
+                              output_grid=output_grid,
+                              method=_method,
+                              atol=_atol,
+                              rtol=_rtol,
+                              first_step=_first_step,
+                              max_step=_max_step)
+
+        sim_phase = SimulationPhase(transcription=tx,
+                                    ode_class=ode_class,
+                                    ode_init_kwargs=ode_init_kwargs)
+
+        sim_phase.time_options.update(self.time_options)
+        sim_phase.time_options['fix_initial'] = True
+        sim_phase.time_options['fix_duration'] = True
+        sim_phase.time_options['initial_bounds'] = (None, None)
+        sim_phase.time_options['duration_bounds'] = (None, None)
+
+        for state_name, state_options in self.state_options.items():
+            sim_phase.state_options[state_name] = deepcopy(state_options)
+            sim_phase.state_options[state_name]['fix_final'] = False
+            sim_phase.state_options[state_name]['input_initial'] = False
+
+        for param_name, param_options in self.parameter_options.items():
+            sim_phase.parameter_options[param_name] = deepcopy(param_options)
+            sim_phase.parameter_options[param_name]['opt'] = False
+
+        for control_name, control_options in self.control_options.items():
+            sim_phase.control_options[control_name] = deepcopy(control_options)
+            sim_phase.control_options[control_name]['opt'] = False
+
+        for pc_name, pc_options in self.polynomial_control_options.items():
+            sim_phase.polynomial_control_options[pc_name] = deepcopy(pc_options)
+            sim_phase.polynomial_control_options[pc_name]['opt'] = False
+
+        sim_phase._timeseries = {ts_name: ts_options for ts_name, ts_options in self._timeseries.items()
+                                 if ts_name == 'timeseries'}
+
+        for state_name, state_options in sim_phase.state_options.items():
+            state_options['fix_final'] = False
+            state_options['input_initial'] = False
+
+        sim_phase.refine_options = deepcopy(self.refine_options)
+        sim_phase.simulate_options = deepcopy(self.simulate_options)
+        sim_phase.timeseries_options = deepcopy(self.timeseries_options)
 
         return sim_phase
 
diff --git a/dymos/phase/simulation_phase.py b/dymos/phase/simulation_phase.py
index 36f33e960..e60fd1c66 100644
--- a/dymos/phase/simulation_phase.py
+++ b/dymos/phase/simulation_phase.py
@@ -17,85 +17,41 @@ class SimulationPhase(Phase):
 
     Parameters
     ----------
-    from_phase : <Phase> or None
-        A phase instance from which the initialized phase should copy its data.
-    times_per_seg : int
-        The number of output points per segment, uniformly distributed.
-    method : str or _unspecified
-        A valid scipy.solve_ivp method for integration.
-    atol : float or _unspecified
-        Absolute error tolerance of the integration.
-    rtol : float or _unspecified
-        Relative error tolerance of the integration.
-    first_step : float or _unspecified
-        Initial step size of the integration.
-    max_step : float or _unspecified
-        Maximum step size of the integration.
-    reports : float or _unspecified
-        If True, generate reports for the subproblem used in integration.
+    transcription : ExplicitShooting
+        The transcription used for the SimulationPhase. It must be an instance of ExplicitShooting.
 
     **kwargs : dict
         Dictionary of optional phase arguments.
     """
-    def __init__(self, from_phase, times_per_seg=_unspecified, method=_unspecified, atol=_unspecified,
-                 rtol=_unspecified, first_step=_unspecified, max_step=_unspecified,
-                 reports=False, **kwargs):
-
-        phase_tx = from_phase.options['transcription']
-        num_seg = phase_tx.grid_data.num_segments
-        seg_order = phase_tx.grid_data.transcription_order
-        seg_ends = phase_tx.grid_data.segment_ends
-        compressed = phase_tx.grid_data.compressed
-
-        sim_options = from_phase.simulate_options
-
-        _method = method if method is not _unspecified else sim_options['method']
-        _atol = atol if atol is not _unspecified else sim_options['atol']
-        _rtol = rtol if rtol is not _unspecified else sim_options['rtol']
-        _first_step = first_step if first_step is not _unspecified else sim_options['first_step']
-        _max_step = max_step if max_step is not _unspecified else sim_options['max_step']
-        _times_per_seg = times_per_seg if times_per_seg is not _unspecified else sim_options['times_per_seg']
-
-        if isinstance(phase_tx, GaussLobatto):
-            grid = GaussLobattoGrid(num_segments=num_seg, nodes_per_seg=seg_order, segment_ends=seg_ends,
-                                    compressed=compressed)
-        elif isinstance(phase_tx, Radau):
-            grid = RadauGrid(num_segments=num_seg, nodes_per_seg=seg_order + 1, segment_ends=seg_ends,
-                             compressed=compressed)
-        elif isinstance(phase_tx.grid_data, GaussLobattoGrid) or \
-                isinstance(phase_tx.grid_data, RadauGrid) or isinstance(phase_tx.grid_data, BirkhoffGrid):
-            grid = phase_tx.grid_data
-        else:
-            raise RuntimeError(f'Unexpected grid class for {phase_tx.grid_data}. Only phases with GaussLobatto '
-                               f'or Radau grids can be simulated.')
-
-        if _times_per_seg is None:
-            output_grid = None
-        else:
-            output_grid = UniformGrid(num_segments=num_seg, nodes_per_seg=_times_per_seg, segment_ends=seg_ends,
-                                      compressed=compressed)
-
-        tx = ExplicitShooting(propagate_derivs=False,
-                              subprob_reports=reports,
-                              grid=grid,
-                              output_grid=output_grid,
-                              method=_method,
-                              atol=_atol,
-                              rtol=_rtol,
-                              first_step=_first_step,
-                              max_step=_max_step)
-
-        super().__init__(from_phase=from_phase, transcription=tx, ode_class=from_phase.options['ode_class'],
-                         ode_init_kwargs=from_phase.options['ode_init_kwargs'])
-
-        # Remove invalid options
-        for state_name, options in self.state_options.items():
-            options['fix_final'] = False  # ExplicitShooting will raise if `fix_final` is True for any states.
-            options['input_initial'] = False  # Only simulate from the initial value, do not connect.
-
-        # Remove all but the default timeseries object
-        self._timeseries = {ts_name: ts_options for ts_name, ts_options in self._timeseries.items()
-                            if ts_name == 'timeseries'}
+    def __init__(self, transcription=None, **kwargs):
+        if not isinstance(transcription, ExplicitShooting):
+            raise ValueError('The transcription for a SimulationPhase must be '
+                             'ExplicitShooting. Use Phase.get_simulation_phase()'
+                             'to create a simulation Phase.')
+        super().__init__(transcription=transcription, **kwargs)
+
+    def duplicate(self, *args, **kwargs):
+        """
+        Create a copy of this phase where most options and attributes are deep copies of those in the original.
+
+        By default, a deepcopy of the transcription in the original phase is used.
+        Boundary constraints, path constraints, and objectives are _NOT_ copied by default, but the user may opt to do so.
+        By default, initial time is not fixed, nor are the initial or final state values.
+        These also can be overridden with the appropriate arguments.
+
+        Parameters
+        ----------
+        *args
+            Additional arguments.
+        **kwargs
+            Keyword arguments.
+
+        Raises
+        ------
+        NotImplmentedError
+            This method is not yet supported for SimulationPhase
+        """
+        raise NotImplementedError('SimulationPhase does not support the duplicate method.')
 
     def set_vals_from_phase(self, from_phase):
         """
@@ -140,72 +96,6 @@ def set_vals_from_phase(self, from_phase):
             val = from_phase.get_val(f'polynomial_controls:{name}', units=options['units'], from_src=False)
             self.set_val(f'polynomial_controls:{name}', val, units=options['units'])
 
-    def initialize_values_from_phase(self, prob, from_phase, phase_path=''):
-        """
-        Initializes values in the Phase using the phase from which it was created.
-
-        Parameters
-        ----------
-        prob : Problem
-            The problem instance to set values taken from the from_phase instance.
-        from_phase : Phase
-            The Phase instance from which the values in this phase are being initialized.
-        phase_path : str
-            The pathname of the system in prob that contains the phases.
-        """
-        phs = from_phase
-
-        op_dict = dict([(name, options) for (name, options) in phs.list_outputs(units=True,
-                                                                                list_autoivcs=True,
-                                                                                out_stream=None)])
-        ip_dict = dict([(name, options) for (name, options) in phs.list_inputs(units=True,
-                                                                               out_stream=None)])
-
-        if self.pathname.partition('.')[0] == self.name:
-            self_path = self.name + '.'
-        else:
-            self_path = self.pathname.partition('.')[0] + '.' + self.name + '.'
-
-        if MPI:
-            op_dict = MPI.COMM_WORLD.bcast(op_dict, root=0)
-
-        # Set the integration times
-        time_name = phs.time_options['name']
-        op = op_dict[f'timeseries.timeseries_comp.{time_name}']
-        prob.set_val(f'{self_path}t_initial', op['val'][0, ...])
-        prob.set_val(f'{self_path}t_duration', op['val'][-1, ...] - op['val'][0, ...])
-
-        # Assign initial state values
-        for name in phs.state_options:
-            prefix = 'states:' if from_phase.timeseries_options['use_prefix'] else ''
-            op = op_dict[f'timeseries.timeseries_comp.{prefix}{name}']
-            prob[f'{self_path}initial_states:{name}'][...] = op['val'][0, ...]
-
-        # Assign control values
-        for name, options in phs.control_options.items():
-            ip = ip_dict[f'control_group.control_interp_comp.controls:{name}']
-            prob[f'{self_path}controls:{name}'][...] = ip['val']
-
-        # Assign polynomial control values
-        for name, options in phs.polynomial_control_options.items():
-            ip = ip_dict[f'polynomial_control_group.interp_comp.'
-                         f'polynomial_controls:{name}']
-            prob[f'{self_path}polynomial_controls:{name}'][...] = ip['val']
-
-        # Assign parameter values
-        for name in phs.parameter_options:
-            units = phs.parameter_options[name]['units']
-
-            # We use this private function to grab the correctly sized variable from the
-            # auto_ivc source.
-            val = phs.get_val(f'parameters:{name}', units=units)
-
-            if phase_path:
-                prob_path = f'{phase_path}.{self.name}.parameters:{name}'
-            else:
-                prob_path = f'{self.name}.parameters:{name}'
-            prob.set_val(prob_path, val)
-
     def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
                                 shape=None, indices=None, lower=None, upper=None, equals=None,
                                 scaler=None, adder=None, ref=None, ref0=None, linear=False, flat_indices=False):
@@ -356,3 +246,55 @@ def add_objective(self, name, loc='final', index=None, shape=(1,), units=None, r
             calculations with other variables sharing the same parallel_deriv_color.
         """
         raise NotImplementedError('SimulationPhase does not support optimization objectives.')
+
+    def check_time_options(self):
+        """
+        Check that time options are valid and issue warnings if invalid options are provided.
+
+        This check is not performed by SimulationPhase.
+
+        Warns
+        -----
+        RuntimeWarning
+            RuntimeWarning is issued in the case of one or more invalid time options.
+        """
+        pass
+
+    def _check_control_options(self):
+        """
+        Check that control options are valid and issue warnings if invalid options are provided.
+
+        This check is not performed by SimulationPhase.
+
+        Warns
+        -----
+        RuntimeWarning
+            RuntimeWarning is issued in the case of one or more invalid time options.
+        """
+        pass
+
+    def _check_polynomial_control_options(self):
+        """
+        Check that polynomial control options are valid and issue warnings if invalid options are provided.
+
+        This check is not performed by SimulationPhase.
+
+        Warns
+        -----
+        RuntimeWarning
+            RuntimeWarning is issued in the case of one or more invalid time options.
+        """
+        pass
+
+    def _check_parameter_options(self):
+        """
+        Check that parameter options are valid and issue warnings if invalid options are provided.
+
+        This check is not performed by SimulationPhase.
+
+        Warns
+        -----
+        RuntimeWarning
+            RuntimeWarning is issued in the case of one or more invalid time options.
+        """
+        pass
diff --git a/dymos/phase/test/test_analytic_phase.py b/dymos/phase/test/test_analytic_phase.py
index 06baafb3d..e47cf6d47 100644
--- a/dymos/phase/test/test_analytic_phase.py
+++ b/dymos/phase/test/test_analytic_phase.py
@@ -150,6 +150,32 @@ def test_simple_bvp(self):
 
         assert_near_equal(y, expected)
 
+    def test_duplicate(self):
+
+        p = om.Problem()
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+        phase = dm.AnalyticPhase(ode_class=SimpleBVPSolution, num_nodes=11)
+
+        phase.set_time_options(units='s', targets=['x'], fix_initial=True, fix_duration=True,
+                               initial_val=0.0, duration_val=1.0)
+        phase.add_parameter('y0', opt=False, val=0.0)
+        phase.add_parameter('y1', opt=False, val=0.0)
+
+        phase2 = phase.duplicate()
+        traj.add_phase('phase2', phase2)
+
+        p.setup()
+
+        p.run_model()
+
+        t = p.get_val('traj.phase2.timeseries.time', units='s')
+        y = p.get_val('traj.phase2.timeseries.y', units='unitless')
+
+        expected = t * (1 - t) * (1 + t - t**2) / 12
+
+        assert_near_equal(y, expected)
+
     def test_renamed_state(self):
 
         class SolutionWithRenamedState(om.ExplicitComponent):
diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py
index 55886fe79..0f4c12db1 100644
--- a/dymos/trajectory/trajectory.py
+++ b/dymos/trajectory/trajectory.py
@@ -60,7 +60,6 @@ def __init__(self, **kwargs):
         self._phase_graph = nx.DiGraph()
         self._has_connected_phases = False
 
-        self.parameter_options = {}
         self.phases = om.ParallelGroup() if self.options['parallel_phases'] else om.Group()
 
     def initialize(self):
@@ -75,6 +74,12 @@ def initialize(self):
                                   '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.')
+        self.options.declare('parameter_options', types=dict, default={},
+                             desc='Options for each parameter in this Trajectory')
+
+    @property
+    def parameter_options(self):
+        return self.options['parameter_options']
 
     def add_phase(self, name, phase, **kwargs):
         """
diff --git a/dymos/visualization/linkage/test/test_gui.py b/dymos/visualization/linkage/test/test_gui.py
index 232b8a8de..641b92472 100644
--- a/dymos/visualization/linkage/test/test_gui.py
+++ b/dymos/visualization/linkage/test/test_gui.py
@@ -9,7 +9,7 @@
 
 if playwright is not None:
     os.system("playwright install")
-    from linkage_report_ui_test import dymos_linkage_gui_test_case  # nopep8: E402
+    from dymos.visualization.linkage.test.linkage_report_ui_test import dymos_linkage_gui_test_case  # nopep8: E402
 
 
 if __name__ == "__main__":
diff --git a/dymos/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py
index 2bdbdc91f..1d67495be 100644
--- a/dymos/visualization/timeseries/bokeh_timeseries_report.py
+++ b/dymos/visualization/timeseries/bokeh_timeseries_report.py
@@ -1,6 +1,7 @@
 import datetime
 from pathlib import Path
 import os.path
+from numpy import isin
 
 try:
     from bokeh.io import output_notebook, output_file, save, show
@@ -16,6 +17,7 @@
 
 import openmdao.api as om
 from openmdao.utils.units import conversion_to_base_units
+from openmdao.recorders.sqlite_recorder import format_version, META_KEY_SEP
 import dymos as dm
 
 
@@ -92,7 +94,7 @@
 """
 
 
-def _meta_tree_subsys_iter(tree, recurse=True, cls=None):
+def _meta_tree_subsys_iter(tree, recurse=True, cls=None, path=None):
     """
     Yield a generator of local subsystems of this system.
 
@@ -102,6 +104,7 @@ def _meta_tree_subsys_iter(tree, recurse=True, cls=None):
         If True, iterate over the whole tree under this system.
     cls : None, str, or Sequence
         The class of the nodes to be iterated
+    path : The absolute path of the given tree.
 
     Yields
     ------
@@ -112,72 +115,179 @@ def _meta_tree_subsys_iter(tree, recurse=True, cls=None):
     for s in tree['children']:
         if s['type'] != 'subsystem':
             continue
+        else:
+            s['path'] = f'{path}.{s["name"]}' if path else s['name']
+
         if cls is None or s['class'] in _cls:
             yield s
         if recurse:
-            for child in _meta_tree_subsys_iter(s, recurse=True, cls=_cls):
+            for child in _meta_tree_subsys_iter(s, recurse=True, cls=_cls, path=s['path']):
                 yield child
 
 
-def _load_data_sources(prob, solution_record_file=None, simulation_record_file=None):
+def _get_model_options(cr, system, run_number=None):
+    """
+    The the options for system stored in the given case reader.
+    If there is more than one set of model options, this function returns the last recorded ones.
+
+    Parameters
+    ----------
+    cr : CaseReader
+        The CaseReader instance holding the data.
+    system : str or None
+        Pathname of system (None for all systems).
+    run_number : int or None
+        Run_driver or run_model iteration to inspect, if given. If None, return the last available.
+
+    Returns
+    -------
+    dict
+        {system: {key: val}}.
+    """
+    if not cr._system_options:
+        raise RuntimeError('System options were not recorded.')
+
+    comp_options = None
+
+    # need to handle edge case for v11 recording
+    if cr._format_version < 12:
+        SEP = '_'
+    else:
+        SEP = META_KEY_SEP
+
+    for key in cr._system_options:
+        if key.find(SEP) > 0:
+            name, num = key.rsplit(SEP, 1)
+        else:
+            name = key
+            num = 0
+
+        # Get the component associated with the highest run number
+        if name == system and (run_number is None or run_number == int(num)):
+            comp_options = cr._system_options[key]['component_options']
+
+    if comp_options is None:
+        raise ValueError(f'No options found for system {system}')
+
+    return comp_options
+
+
+def _get_trajs_and_phases(cr):
+    trajs = {}
+    traj_cls = 'dymos.trajectory.trajectory:Trajectory'
+    phase_cls = ['dymos.phase.phase:Phase',
+                 'dymos.phase.phase:AnalyticPhase',
+                 'dymos.phase.phase:SimualationPhase']
+
+    traj_nodes = [n for n in _meta_tree_subsys_iter(cr.problem_metadata['tree'], cls=traj_cls)]
+    phase_nodes = [n for n in _meta_tree_subsys_iter(cr.problem_metadata['tree'], cls=phase_cls)]
+
+    for tn in traj_nodes:
+        phase_nodes = [n for n in _meta_tree_subsys_iter(tn, cls=phase_cls, path=tn['path'])]
+        traj_path = tn['path']
+        trajs[traj_path] = {'phases': {}, 'parameter_options': {}, 'name': tn['name']}
+        traj_options = _get_model_options(cr=cr, system=traj_path)
+        for param_name, param_options in traj_options['parameter_options'].items():
+            trajs[traj_path]['parameter_options'][param_name] = param_options
+        for pn in phase_nodes:
+            phase_path = pn['path']
+            phase_options = _get_model_options(cr=cr, system=phase_path)
+            phase_meta = trajs[traj_path]['phases'][phase_path] = {'time_options': None,
+                                                                   'parameter_options': None,
+                                                                   'state_options': None,
+                                                                   'control_options': None,
+                                                                   'polynomial_control_options': None,
+                                                                   'name': pn['name']}
+            phase_meta['time_options'] = phase_options['time_options']
+            phase_meta['parameter_options'] = phase_options['parameter_options']
+            phase_meta['state_options'] = phase_options['state_options']
+            phase_meta['control_options'] = phase_options['control_options']
+            phase_meta['polynomial_control_options'] = phase_options['polynomial_control_options']
+
+    return trajs
+
+
+def _load_data_sources(solution_record_file=None, simulation_record_file=None):
+    """
+    Load the data for the timeseries plots from the given solution and record files.
+
+    Parameters
+    ----------
+    solution_record_file : str
+        The path to the solution record file.
+    sim_record_file : str
+        The path to the corresponding simulation record file.
+
+    Returns
+    -------
+    dict
+        A dictionary containing parameter data, solution timeseries data, simulation timeseries data, and
+        units for each timeseries output.
+    """
 
+    traj_and_phase_meta = None
     data_dict = {}
 
     if Path(solution_record_file).is_file():
         sol_cr = om.CaseReader(solution_record_file)
         sol_case = sol_cr.get_case('final')
-        outputs = {name: meta for name, meta in sol_case.list_outputs(units=True, out_stream=None)}
         abs2prom_map = sol_cr.problem_metadata['abs2prom']
+        traj_and_phase_meta = _get_trajs_and_phases(sol_cr)
     else:
         sol_case = None
 
     if Path(simulation_record_file).is_file():
         sim_cr = om.CaseReader(simulation_record_file)
         sim_case = sim_cr.get_case('final')
-        outputs = {name: meta for name, meta in sim_case.list_outputs(units=True, out_stream=None)}
         abs2prom_map = sim_cr.problem_metadata['abs2prom']
+        if traj_and_phase_meta is None:
+            traj_and_phase_meta = _get_trajs_and_phases(sim_cr)
     else:
         sim_case = None
 
+    source_case = sol_case if sol_case else sim_case
+    outputs = {abs_path: meta for abs_path, meta in source_case.list_outputs(out_stream=None, units=True)}
+
     if sol_cr is None and sim_cr is None:
         om.issue_warning('No recorded data provided. Trajectory results report will not be created.')
         return
 
-    for traj in prob.model.system_iter(include_self=True, recurse=True, typ=dm.Trajectory):
-        traj_name = traj.pathname.split('.')[-1]
-        data_dict[traj_name] = {'param_data_by_phase': {},
-                                'sol_data_by_phase': {},
-                                'sim_data_by_phase': {},
-                                'timeseries_units': {}}
+    for traj_path, traj_data in traj_and_phase_meta.items():
+        traj_params = traj_data['parameter_options']
+        traj_name = traj_data['name']
+        data_dict[traj_data['name']] = {'param_data_by_phase': {},
+                                        'sol_data_by_phase': {},
+                                        'sim_data_by_phase': {},
+                                        'timeseries_units': {}}
 
-        for phase_name, phase in traj._phases.items():
+        for phase_path, phase_data in traj_data['phases'].items():
+            phase_name = phase_data['name']
+            phase_param_data = \
+                data_dict[traj_name]['param_data_by_phase'][phase_name] = \
+                {'param': [], 'val': [], 'units': []}
 
-            data_dict[traj_name]['param_data_by_phase'][phase_name] = {'param': [], 'val': [], 'units': []}
-            phase_sol_data = data_dict[traj_name]['sol_data_by_phase'][phase_name] = {}
-            phase_sim_data = data_dict[traj_name]['sim_data_by_phase'][phase_name] = {}
-            ts_units_dict = data_dict[traj_name]['timeseries_units']
+            phase_sol_data = data_dict[traj_data['name']]['sol_data_by_phase'][phase_name] = {}
+            phase_sim_data = data_dict[traj_data['name']]['sim_data_by_phase'][phase_name] = {}
 
             param_outputs = {op: meta for op, meta in outputs.items()
-                             if op.startswith(f'{phase.pathname}.param_comp.parameter_vals')}
-            param_case = sol_case if sol_case else sim_case
-
-            for output_name in sorted(param_outputs.keys(), key=str.casefold):
-                meta = param_outputs[output_name]
-                param_dict = data_dict[traj_name]['param_data_by_phase'][phase_name]
-
-                prom_name = abs2prom_map['output'][output_name]
-                param_name = output_name.replace(f'{phase.pathname}.param_comp.parameter_vals:', '', 1)
-
-                param_dict['param'].append(param_name)
-                param_dict['units'].append(meta['units'])
-                param_dict['val'].append(param_case.get_val(prom_name, units=meta['units']))
-
-            ts_outputs = {op: meta for op, meta in outputs.items() if op.startswith(f'{phase.pathname}.timeseries')}
+                             if op.startswith(f'{phase_path}.param_comp.parameter_vals:')}
+            ts_outputs = {op: meta for op, meta in outputs.items()
+                          if op.startswith(f'{phase_path}.timeseries.')}
+
+            # Populate the phase parameter data
+            phase_params = traj_and_phase_meta[traj_path]['phases'][phase_path]['parameter_options']
+            for param_name in sorted(phase_params, key=str.casefold):
+                units = phase_params[param_name]['units']
+                param_path = f'{traj_path}.{phase_name}.parameter_vals:{param_name}'
+                phase_param_data['param'].append(param_name)
+                phase_param_data['units'].append(units)
+                phase_param_data['val'].append(source_case.get_val(param_path, units=units))
 
             # Find the "largest" unit used for any timeseries output across all phases
-            for output_name in sorted(ts_outputs.keys(), key=str.casefold):
-                meta = ts_outputs[output_name]
-                prom_name = abs2prom_map['output'][output_name]
+            ts_units_dict = data_dict[traj_data['name']]['timeseries_units']
+            for abs_name in sorted(ts_outputs.keys(), key=str.casefold):
+                meta = ts_outputs[abs_name]
+                prom_name = meta['prom_name']
                 var_name = prom_name.split('.')[-1]
 
                 if var_name not in ts_units_dict:
@@ -188,16 +298,7 @@ def _load_data_sources(prob, solution_record_file=None, simulation_record_file=N
                     if new_conv_factor < old_conv_factor:
                         ts_units_dict[var_name] = meta['units']
 
-        # Now a second pass through the phases since we know the units in which to plot
-        # each timeseries variable output.
-        for phase_name, phase in traj._phases.items():
-
-            phase_sol_data = data_dict[traj_name]['sol_data_by_phase'][phase_name] = {}
-            phase_sim_data = data_dict[traj_name]['sim_data_by_phase'][phase_name] = {}
-            ts_units_dict = data_dict[traj_name]['timeseries_units']
-
-            ts_outputs = {op: meta for op, meta in outputs.items() if op.startswith(f'{phase.pathname}.timeseries')}
-
+            # Populate the phase timeseries data
             for output_name in sorted(ts_outputs.keys(), key=str.casefold):
                 meta = ts_outputs[output_name]
                 prom_name = abs2prom_map['output'][output_name]
@@ -233,23 +334,23 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi
     theme : str
         A valid bokeh theme name to style the report.
     """
+    report_dir = Path(prob.get_reports_dir()) if prob is not None else Path(os.getcwd())
+    if not report_dir.exists():
+        report_dir = Path(os.getcwd())
+
     # For the primary timeseries in each phase in each trajectory, build a set of the pathnames
     # to be plotted.
-    source_data = _load_data_sources(prob, solution_record_file, simulation_record_file)
+    source_data = _load_data_sources(solution_record_file, simulation_record_file)
 
     # Colors of each phase in the plot. Start with the bright colors followed by the faded ones.
     if not _NO_BOKEH:
         colors = bp.d3['Category20'][20][0::2] + bp.d3['Category20'][20][1::2]
         curdoc().theme = theme
 
-    for traj in prob.model.system_iter(include_self=True, recurse=True, typ=dm.Trajectory):
-        traj_name = traj.pathname.split('.')[-1]
-        report_filename = f'{traj.pathname}_results_report.html'
-        report_dir = Path(prob.get_reports_dir())
+    for traj_path, traj_data in source_data.items():
+        traj_name = traj_path.split('.')[-1]
+        report_filename = f'{traj_path}_results_report.html'
         report_path = report_dir / report_filename
-        if not os.path.isdir(report_dir):
-            om.issue_warning(f'Reports directory not available. {report_path} will not be created.')
-            continue
         if _NO_BOKEH:
             with open(report_path, 'wb') as f:
                 f.write("<html>\n<head>\n<title> \nError: bokeh not available</title>\n</head> <body>\n"
@@ -258,21 +359,19 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi
             continue
 
         param_tables = []
-        phase_names = []
-
-        for phase_name, phase in traj._phases.items():
-
-            phase_names.append(phase_name)
-
+        phase_names = [k.split('.')[-1] for k in traj_data['param_data_by_phase']]
+        for phase_name, param_data in traj_data['param_data_by_phase'].items():
             # Make the parameter table
-            source = ColumnDataSource(source_data[traj_name]['param_data_by_phase'][phase_name])
             columns = [
                 TableColumn(field='param', title='Parameter'),
                 TableColumn(field='val', title='Value'),
                 TableColumn(field='units', title='Units'),
             ]
-            param_tables.append(DataTable(source=source, columns=columns, index_position=None,
-                                          height=30*len(source.data['param']), sizing_mode='stretch_both'))
+            param_tables.append(DataTable(source=ColumnDataSource(param_data),
+                                          columns=columns,
+                                          index_position=None,
+                                          height=30*len(param_data['param']),
+                                          sizing_mode='stretch_both'))
 
         # Plot the timeseries
         ts_units_dict = source_data[traj_name]['timeseries_units']
diff --git a/joss/flow_charts/coupled_co_design.tex b/joss/flow_charts/coupled_co_design.tex
new file mode 100644
index 000000000..3f88ebd56
--- /dev/null
+++ b/joss/flow_charts/coupled_co_design.tex
@@ -0,0 +1,27 @@
+
+% XDSM diagram created with pyXDSM 2.1.3.
+\documentclass{article}
+\usepackage{geometry}
+\usepackage{amsfonts}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{tikz}
+
+% Optional packages such as sfmath set through python interface
+\usepackage{}
+
+% Define the set of TikZ packages to be included in the architecture diagram document
+\usetikzlibrary{arrows,chains,positioning,scopes,shapes.geometric,shapes.misc,shadows}
+
+
+% Set the border around all of the architecture diagrams to be tight to the diagrams themselves
+% (i.e. no longer need to tinker with page size parameters)
+\usepackage[active,tightpage]{preview}
+\PreviewEnvironment{tikzpicture}
+\setlength{\PreviewBorder}{5pt}
+
+\begin{document}
+
+\input{"coupled_co_design.tikz"}
+
+\end{document}
diff --git a/joss/flow_charts/coupled_co_design.tikz b/joss/flow_charts/coupled_co_design.tikz
new file mode 100644
index 000000000..891a136e4
--- /dev/null
+++ b/joss/flow_charts/coupled_co_design.tikz
@@ -0,0 +1,55 @@
+
+%%% Preamble Requirements %%%
+% \usepackage{geometry}
+% \usepackage{amsfonts}
+% \usepackage{amsmath}
+% \usepackage{amssymb}
+% \usepackage{tikz}
+
+% Optional packages such as sfmath set through python interface
+% \usepackage{}
+
+% \usetikzlibrary{arrows,chains,positioning,scopes,shapes.geometric,shapes.misc,shadows}
+
+%%% End Preamble Requirements %%%
+
+\input{"/Users/rfalck/anaconda3/envs/py37/lib/python3.7/site-packages/pyxdsm/diagram_styles"}
+\begin{tikzpicture}
+
+\matrix[MatrixSetup]{
+%Row 0
+\node [Optimization] (OPT) {$\text{Optimizer}$};&
+\node [DataInter] (OPT-static) {$\begin{array}{c}\bar{d}\end{array}$};&
+\node [DataInter] (OPT-dynamic) {$\begin{array}{c}t, \bar{x}, \bar{u}\end{array}$};\\
+%Row 1
+\node [DataInter] (static-OPT) {$\begin{array}{c}g_\text{static}\end{array}$};&
+\node [Group] (static) {$\text{Static System Model}$};&
+\node [DataInter] (static-dynamic) {$\begin{array}{c}\text{static} \\ \text{outputs}\end{array}$};\\
+%Row 2
+\node [DataInter] (dynamic-OPT) {$\begin{array}{c}J, \bar{g}_0, \bar{g}_f, \bar{p}\end{array}$};&
+\node [DataInter] (dynamic-static) {$\begin{array}{c}\text{dynamic} \\ \text{outputs}\end{array}$};&
+\node [Group] (dynamic) {$\text{ODE or DAE}$};\\
+};
+
+% XDSM process chains
+
+
+\begin{pgfonlayer}{data}
+\path
+% Horizontal edges
+(dynamic) edge [DataLine] (dynamic-OPT)
+(OPT) edge [DataLine] (OPT-dynamic)
+(static) edge [DataLine] (static-OPT)
+(OPT) edge [DataLine] (OPT-static)
+(static) edge [DataLine] (static-dynamic)
+(dynamic) edge [DataLine] (dynamic-static)
+% Vertical edges
+(dynamic-OPT) edge [DataLine] (OPT)
+(OPT-dynamic) edge [DataLine] (dynamic)
+(static-OPT) edge [DataLine] (OPT)
+(OPT-static) edge [DataLine] (static)
+(static-dynamic) edge [DataLine] (dynamic)
+(dynamic-static) edge [DataLine] (static);
+\end{pgfonlayer}
+
+\end{tikzpicture}
diff --git a/joss/flow_charts/opt_control.tex b/joss/flow_charts/opt_control.tex
new file mode 100644
index 000000000..aec07cbcf
--- /dev/null
+++ b/joss/flow_charts/opt_control.tex
@@ -0,0 +1,27 @@
+
+% XDSM diagram created with pyXDSM 2.1.3.
+\documentclass{article}
+\usepackage{geometry}
+\usepackage{amsfonts}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{tikz}
+
+% Optional packages such as sfmath set through python interface
+\usepackage{}
+
+% Define the set of TikZ packages to be included in the architecture diagram document
+\usetikzlibrary{arrows,chains,positioning,scopes,shapes.geometric,shapes.misc,shadows}
+
+
+% Set the border around all of the architecture diagrams to be tight to the diagrams themselves
+% (i.e. no longer need to tinker with page size parameters)
+\usepackage[active,tightpage]{preview}
+\PreviewEnvironment{tikzpicture}
+\setlength{\PreviewBorder}{5pt}
+
+\begin{document}
+
+\input{"opt_control.tikz"}
+
+\end{document}
diff --git a/joss/flow_charts/opt_control.tikz b/joss/flow_charts/opt_control.tikz
new file mode 100644
index 000000000..2c6e777b2
--- /dev/null
+++ b/joss/flow_charts/opt_control.tikz
@@ -0,0 +1,41 @@
+
+%%% Preamble Requirements %%%
+% \usepackage{geometry}
+% \usepackage{amsfonts}
+% \usepackage{amsmath}
+% \usepackage{amssymb}
+% \usepackage{tikz}
+
+% Optional packages such as sfmath set through python interface
+% \usepackage{}
+
+% \usetikzlibrary{arrows,chains,positioning,scopes,shapes.geometric,shapes.misc,shadows}
+
+%%% End Preamble Requirements %%%
+
+\input{"/Users/rfalck/anaconda3/envs/py37/lib/python3.7/site-packages/pyxdsm/diagram_styles"}
+\begin{tikzpicture}
+
+\matrix[MatrixSetup]{
+%Row 0
+\node [Optimization] (OPT) {$\text{Optimizer}$};&
+\node [DataInter] (OPT-ODE) {$\begin{array}{c}t, \bar{x}, \bar{u}, \bar{d}\end{array}$};\\
+%Row 1
+\node [DataInter] (ODE-OPT) {$\begin{array}{c}J, \bar{g}_0, \bar{g}_f, \bar{p}\end{array}$};&
+\node [Group] (ODE) {$\text{ODE or DAE}$};\\
+};
+
+% XDSM process chains
+
+
+\begin{pgfonlayer}{data}
+\path
+% Horizontal edges
+(ODE) edge [DataLine] (ODE-OPT)
+(OPT) edge [DataLine] (OPT-ODE)
+% Vertical edges
+(ODE-OPT) edge [DataLine] (OPT)
+(OPT-ODE) edge [DataLine] (ODE);
+\end{pgfonlayer}
+
+\end{tikzpicture}
diff --git a/joss/flow_charts/sequential_co_design.tex b/joss/flow_charts/sequential_co_design.tex
new file mode 100644
index 000000000..b7696981f
--- /dev/null
+++ b/joss/flow_charts/sequential_co_design.tex
@@ -0,0 +1,27 @@
+
+% XDSM diagram created with pyXDSM 2.1.3.
+\documentclass{article}
+\usepackage{geometry}
+\usepackage{amsfonts}
+\usepackage{amsmath}
+\usepackage{amssymb}
+\usepackage{tikz}
+
+% Optional packages such as sfmath set through python interface
+\usepackage{}
+
+% Define the set of TikZ packages to be included in the architecture diagram document
+\usetikzlibrary{arrows,chains,positioning,scopes,shapes.geometric,shapes.misc,shadows}
+
+
+% Set the border around all of the architecture diagrams to be tight to the diagrams themselves
+% (i.e. no longer need to tinker with page size parameters)
+\usepackage[active,tightpage]{preview}
+\PreviewEnvironment{tikzpicture}
+\setlength{\PreviewBorder}{5pt}
+
+\begin{document}
+
+\input{"sequential_co_design.tikz"}
+
+\end{document}
diff --git a/joss/flow_charts/sequential_co_design.tikz b/joss/flow_charts/sequential_co_design.tikz
new file mode 100644
index 000000000..b8808f526
--- /dev/null
+++ b/joss/flow_charts/sequential_co_design.tikz
@@ -0,0 +1,63 @@
+
+%%% Preamble Requirements %%%
+% \usepackage{geometry}
+% \usepackage{amsfonts}
+% \usepackage{amsmath}
+% \usepackage{amssymb}
+% \usepackage{tikz}
+
+% Optional packages such as sfmath set through python interface
+% \usepackage{}
+
+% \usetikzlibrary{arrows,chains,positioning,scopes,shapes.geometric,shapes.misc,shadows}
+
+%%% End Preamble Requirements %%%
+
+\input{"/Users/rfalck/anaconda3/envs/py37/lib/python3.7/site-packages/pyxdsm/diagram_styles"}
+\begin{tikzpicture}
+
+\matrix[MatrixSetup]{
+%Row 0
+\node [Optimization] (OPT_static) {$\text{Static Optimizer}$};&
+\node [DataInter] (OPT_static-static) {$\begin{array}{c}\bar{d}\end{array}$};&
+\node [DataInter] (OPT_static-OPT_dynamic) {$\begin{array}{c}\text{static optimization} \\ \text{outputs}\end{array}$};&
+\\
+%Row 1
+\node [DataInter] (static-OPT_static) {$\begin{array}{c}J_\text{static}, g_\text{static}\end{array}$};&
+\node [Group] (static) {$\text{Static System Model}$};&
+&
+\\
+%Row 2
+\node [DataInter] (OPT_dynamic-OPT_static) {$\begin{array}{c}\text{dynamic optimization} \\ \text{outputs}\end{array}$};&
+&
+\node [Optimization] (OPT_dynamic) {$\text{Dynamic Optimizer}$};&
+\node [DataInter] (OPT_dynamic-dynamic) {$\begin{array}{c}t, \bar{x}, \bar{u}\end{array}$};\\
+%Row 3
+&
+&
+\node [DataInter] (dynamic-OPT_dynamic) {$\begin{array}{c}J_\text{dynamic}, \bar{g}_0, \bar{g}_f, \bar{p}\end{array}$};&
+\node [Group] (dynamic) {$\text{ODE or DAE}$};\\
+};
+
+% XDSM process chains
+
+
+\begin{pgfonlayer}{data}
+\path
+% Horizontal edges
+(dynamic) edge [DataLine] (dynamic-OPT_dynamic)
+(OPT_dynamic) edge [DataLine] (OPT_dynamic-dynamic)
+(static) edge [DataLine] (static-OPT_static)
+(OPT_static) edge [DataLine] (OPT_static-static)
+(OPT_static) edge [DataLine] (OPT_static-OPT_dynamic)
+(OPT_dynamic) edge [DataLine] (OPT_dynamic-OPT_static)
+% Vertical edges
+(dynamic-OPT_dynamic) edge [DataLine] (OPT_dynamic)
+(OPT_dynamic-dynamic) edge [DataLine] (dynamic)
+(static-OPT_static) edge [DataLine] (OPT_static)
+(OPT_static-static) edge [DataLine] (static)
+(OPT_static-OPT_dynamic) edge [DataLine] (OPT_dynamic)
+(OPT_dynamic-OPT_static) edge [DataLine] (OPT_static);
+\end{pgfonlayer}
+
+\end{tikzpicture}
diff --git a/joss/test/test_cannonball_for_joss.py b/joss/test/test_cannonball_for_joss.py
index 472fde06d..395521135 100644
--- a/joss/test/test_cannonball_for_joss.py
+++ b/joss/test/test_cannonball_for_joss.py
@@ -3,7 +3,7 @@
 from openmdao.utils.assert_utils import assert_near_equal
 
 
-@use_tempdirs
+# @use_tempdirs
 class TestCannonballForJOSS(unittest.TestCase):
 
     @require_pyoptsparse(optimizer='SLSQP')
@@ -11,6 +11,7 @@ def test_results(self):
         # begin code for paper
         import numpy as np
         from scipy.interpolate import interp1d
+        import matplotlib
         import matplotlib.pyplot as plt
 
         import openmdao.api as om
@@ -222,7 +223,7 @@ def compute(self, inputs, outputs):
         y0 = p.get_val('traj.ascent.timeseries.h', units='m')
         x1 = p.get_val('traj.descent.timeseries.r', units='m')
         y1 = p.get_val('traj.descent.timeseries.h', units='m')
-        tab20 = plt.cm.get_cmap('tab20').colors
+        tab20 = matplotlib.colormaps['tab20'].colors
         ax.plot(x0, y0, marker='o', label='ascent', color=tab20[0])
         ax.plot(x1, y1, marker='o', label='descent', color=tab20[1])
         ax.legend(loc='best')
@@ -232,3 +233,7 @@ def compute(self, inputs, outputs):
         # End code for paper
 
         assert_near_equal(x1[-1], 3064, tolerance=1.0E-4)
+
+
+if __name__ == '__main__':
+    unittest.main()