diff --git a/docs/dymos_book/_toc.yml b/docs/dymos_book/_toc.yml
index 29561f73f..9d3147daf 100644
--- a/docs/dymos_book/_toc.yml
+++ b/docs/dymos_book/_toc.yml
@@ -41,6 +41,7 @@ parts:
     - file: examples/water_rocket/water_rocket
     - file: examples/cart_pole/cart_pole
     - file: examples/brachistochrone/brachistochrone_tandem_phases
+    - file: examples/cannonball_implicit_duration/cannonball_implicit_duration
 - caption: Feature Reference
   chapters:
   - file: features/phases/phases_list
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
new file mode 100644
index 000000000..2eab4113c
--- /dev/null
+++ b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb
@@ -0,0 +1,471 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# This cell is mandatory in all Dymos documentation notebooks.\n",
+    "missing_packages = []\n",
+    "try:\n",
+    "    import openmdao.api as om\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !python -m pip install openmdao[notebooks]\n",
+    "    else:\n",
+    "        missing_packages.append('openmdao')\n",
+    "try:\n",
+    "    import dymos as dm\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !python -m pip install dymos\n",
+    "    else:\n",
+    "        missing_packages.append('dymos')\n",
+    "try:\n",
+    "    import pyoptsparse\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !pip install -q condacolab\n",
+    "        import condacolab\n",
+    "        condacolab.install_miniconda()\n",
+    "        !conda install -c conda-forge pyoptsparse\n",
+    "    else:\n",
+    "        missing_packages.append('pyoptsparse')\n",
+    "if missing_packages:\n",
+    "    raise EnvironmentError('This notebook requires the following packages '\n",
+    "                           'please install them and restart this notebook\\'s runtime: {\",\".join(missing_packages)}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Cannonball Implicit Duration\n",
+    "\n",
+    "This example demonstrates determining the duration of a phase using a \n",
+    "nonlinear solver to satisfy an endpoint condition. As in the multi-phase\n",
+    "cannonball problem, we will simultaneously find the optimal design for \n",
+    "the cannonball (its radius) and the optimal flight profile (its launch\n",
+    "angle). However, here we will use a single phase that terminates at the\n",
+    "condition\n",
+    "\n",
+    "\\begin{align}\n",
+    "    h(t_f) = 0\n",
+    "\\end{align}\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The cannonball sizing component\n",
+    "\n",
+    "This component computes the area (needed to compute drag) and mass (needed to compute energy) of a cannonball of a given radius and density.\n",
+    "\n",
+    "This component sits upstream of the trajectory model and feeds its outputs to the trajectory as parameters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "output_scroll"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from scipy.interpolate import interp1d\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "import openmdao.api as om\n",
+    "\n",
+    "import dymos as dm\n",
+    "from dymos.models.atmosphere.atmos_1976 import USatm1976Data\n",
+    "\n",
+    "#############################################\n",
+    "# Component for the design part of the model\n",
+    "#############################################\n",
+    "class CannonballSizeComp(om.ExplicitComponent):\n",
+    "    \"\"\"\n",
+    "    Compute the area and mass of a cannonball with a given radius and density.\n",
+    "\n",
+    "    Notes\n",
+    "    -----\n",
+    "    This component is not vectorized with 'num_nodes' as is the usual way\n",
+    "    with Dymos, but is instead intended to compute a scalar mass and reference\n",
+    "    area from scalar radius and density inputs. This component does not reside\n",
+    "    in the ODE but instead its outputs are connected to the trajectory via\n",
+    "    input design parameters.\n",
+    "    \"\"\"\n",
+    "    def setup(self):\n",
+    "        self.add_input(name='radius', val=1.0, desc='cannonball radius', units='m')\n",
+    "        self.add_input(name='dens', val=7870., desc='cannonball density', units='kg/m**3')\n",
+    "\n",
+    "        self.add_output(name='mass', shape=(1,), desc='cannonball mass', units='kg')\n",
+    "        self.add_output(name='S', shape=(1,), desc='aerodynamic reference area', units='m**2')\n",
+    "\n",
+    "        self.declare_partials(of='mass', wrt='dens')\n",
+    "        self.declare_partials(of='mass', wrt='radius')\n",
+    "\n",
+    "        self.declare_partials(of='S', wrt='radius')\n",
+    "\n",
+    "    def compute(self, inputs, outputs):\n",
+    "        radius = inputs['radius']\n",
+    "        dens = inputs['dens']\n",
+    "\n",
+    "        outputs['mass'] = (4/3.) * dens * np.pi * radius ** 3\n",
+    "        outputs['S'] = np.pi * radius ** 2\n",
+    "\n",
+    "    def compute_partials(self, inputs, partials):\n",
+    "        radius = inputs['radius']\n",
+    "        dens = inputs['dens']\n",
+    "\n",
+    "        partials['mass', 'dens'] = (4/3.) * np.pi * radius ** 3\n",
+    "        partials['mass', 'radius'] = 4. * dens * np.pi * radius ** 2\n",
+    "\n",
+    "        partials['S', 'radius'] = 2 * np.pi * radius\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The cannonball ODE component\n",
+    "\n",
+    "This component computes the state rates and the kinetic energy of the cannonball.\n",
+    "By calling the `declare_coloring` method wrt all inputs and using method `'cs'`, we're telling OpenMDAO to automatically determine the sparsity pattern of the outputs with respect to the inputs, **and** to automatically compute those outputs using complex-step approximation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class CannonballODE(om.ExplicitComponent):\n",
+    "    \"\"\"\n",
+    "    Cannonball ODE assuming flat earth and accounting for air resistance\n",
+    "    \"\"\"\n",
+    "\n",
+    "    def initialize(self):\n",
+    "        self.options.declare('num_nodes', types=int)\n",
+    "\n",
+    "    def setup(self):\n",
+    "        nn = self.options['num_nodes']\n",
+    "\n",
+    "        # static parameters\n",
+    "        self.add_input('m', units='kg')\n",
+    "        self.add_input('S', units='m**2')\n",
+    "        # 0.5 good assumption for a sphere\n",
+    "        self.add_input('CD', 0.5)\n",
+    "\n",
+    "        # time varying inputs\n",
+    "        self.add_input('h', units='m', shape=nn)\n",
+    "        self.add_input('v', units='m/s', shape=nn)\n",
+    "        self.add_input('gam', units='rad', shape=nn)\n",
+    "\n",
+    "        # state rates\n",
+    "        self.add_output('v_dot', shape=nn, units='m/s**2', tags=['dymos.state_rate_source:v'])\n",
+    "        self.add_output('gam_dot', shape=nn, units='rad/s', tags=['dymos.state_rate_source:gam'])\n",
+    "        self.add_output('h_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:h'])\n",
+    "        self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r'])\n",
+    "        self.add_output('ke', shape=nn, units='J')\n",
+    "\n",
+    "        # Ask OpenMDAO to compute the partial derivatives using complex-step\n",
+    "        # with a partial coloring algorithm for improved performance, and use\n",
+    "        # a graph coloring algorithm to automatically detect the sparsity pattern.\n",
+    "        self.declare_coloring(wrt='*', method='cs')\n",
+    "\n",
+    "        alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]\n",
+    "        rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]\n",
+    "        self.rho_interp = interp1d(np.array(alt_data, dtype=complex),\n",
+    "                                   np.array(rho_data, dtype=complex),\n",
+    "                                   kind='linear')\n",
+    "\n",
+    "    def compute(self, inputs, outputs):\n",
+    "\n",
+    "        gam = inputs['gam']\n",
+    "        v = inputs['v']\n",
+    "        h = inputs['h']\n",
+    "        m = inputs['m']\n",
+    "        S = inputs['S']\n",
+    "        CD = inputs['CD']\n",
+    "\n",
+    "        GRAVITY = 9.80665  # m/s**2\n",
+    "\n",
+    "        # handle complex-step gracefully from the interpolant\n",
+    "        if np.iscomplexobj(h):\n",
+    "            rho = self.rho_interp(inputs['h'])\n",
+    "        else:\n",
+    "            rho = self.rho_interp(inputs['h']).real\n",
+    "\n",
+    "        q = 0.5*rho*inputs['v']**2\n",
+    "        qS = q * S\n",
+    "        D = qS * CD\n",
+    "        cgam = np.cos(gam)\n",
+    "        sgam = np.sin(gam)\n",
+    "        outputs['v_dot'] = - D/m-GRAVITY*sgam\n",
+    "        outputs['gam_dot'] = -(GRAVITY/v)*cgam\n",
+    "        outputs['h_dot'] = v*sgam\n",
+    "        outputs['r_dot'] = v*cgam\n",
+    "        outputs['ke'] = 0.5*m*v**2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generating an initial guess\n",
+    "\n",
+    "Since we will be using a nonlinear solver to converge the dynamics and the duration, a simple linear initial guess will not be sufficient. The initial guess we will use is generated from the kinematics equations that assume no atmospheric drag. We compute the state values for some duration and initial flight path angle\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def initial_guess(t_dur, gam0=np.pi/3):\n",
+    "    t = np.linspace(0, t_dur, num=100)\n",
+    "    g = -9.80665\n",
+    "    v0 = -0.5*g*t_dur/np.sin(gam0)\n",
+    "    r = v0*np.cos(gam0)*t\n",
+    "    h = v0*np.sin(gam0)*t + 0.5*g*t**2\n",
+    "    v = np.sqrt(v0*np.cos(gam0)**2 + (v0*np.sin(gam0) + g*t)**2)\n",
+    "    gam = np.arctan2(v0*np.sin(gam0) + g*t, v0*np.cos(gam0)**2)\n",
+    "\n",
+    "    guess = {'t': t, 'r': r, 'h': h, 'v': v, 'gam': gam}\n",
+    "\n",
+    "    return guess\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Building and running the problem\n",
+    "\n",
+    "The following code defines the components for the physical\n",
+    "cannonball calculations and ODE problem, sets up trajectory using a\n",
+    "single phase, and specifies the stopping condition for the phase.\n",
+    "The initial flight path angle is free, since 45 degrees is not\n",
+    "necessarily optimal once air resistance is taken into account."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p = om.Problem(model=om.Group())\n",
+    "\n",
+    "p.driver = om.pyOptSparseDriver()\n",
+    "p.driver.options['optimizer'] = 'IPOPT'\n",
+    "p.driver.declare_coloring()\n",
+    "\n",
+    "p.driver.opt_settings['derivative_test'] = 'first-order'\n",
+    "p.driver.opt_settings['mu_strategy'] = 'monotone'\n",
+    "p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'\n",
+    "p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'\n",
+    "p.driver.opt_settings['mu_init'] = 0.01\n",
+    "p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'\n",
+    "\n",
+    "p.set_solver_print(level=0, depth=99)\n",
+    "\n",
+    "p.model.add_subsystem('size_comp', CannonballSizeComp(),\n",
+    "                      promotes_inputs=['radius', 'dens'])\n",
+    "p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')\n",
+    "p.model.add_design_var('radius', lower=0.01, upper=0.10,\n",
+    "                       ref0=0.01, ref=0.10, units='m')\n",
+    "\n",
+    "traj = p.model.add_subsystem('traj', dm.Trajectory())\n",
+    "\n",
+    "transcription = dm.Radau(num_segments=25, order=3, compressed=False)\n",
+    "phase = dm.Phase(ode_class=CannonballODE, transcription=transcription)\n",
+    "\n",
+    "phase = traj.add_phase('phase', phase)\n",
+    "\n",
+    "# Duration is bounded to be greater than 1 to ensure the nonlinear solve\n",
+    "# does not converge to the initial point.\n",
+    "phase.set_time_options(fix_initial=True, duration_bounds=(1, 100), units='s')\n",
+    "\n",
+    "# All initial states except flight path angle are fixed\n",
+    "# The output of the ODE which provides the rate source for each state\n",
+    "# is obtained from the tags used on those outputs in the ODE.\n",
+    "# The units of the states are automatically inferred by multiplying the units\n",
+    "# of those rates by the time units.\n",
+    "phase.add_state('r', fix_initial=True, solve_segments='forward')\n",
+    "phase.add_state('h', fix_initial=True, solve_segments='forward')\n",
+    "phase.add_state('gam', fix_initial=False, solve_segments='forward')\n",
+    "phase.add_state('v', fix_initial=False, solve_segments='forward')\n",
+    "\n",
+    "phase.add_parameter('S', units='m**2', static_target=True)\n",
+    "phase.add_parameter('m', units='kg', static_target=True)\n",
+    "phase.add_parameter('CD', val=0.5, opt=False, static_target=True)\n",
+    "\n",
+    "phase.add_boundary_constraint('ke', loc='initial',\n",
+    "                              upper=400000, lower=0, ref=100000)\n",
+    "\n",
+    "# A duration balance is added setting altitude to zero.\n",
+    "# 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",
+    "phase.add_objective('r', loc='final', scaler=-1.0)\n",
+    "\n",
+    "p.model.connect('size_comp.mass', 'traj.phase.parameters:m')\n",
+    "p.model.connect('size_comp.S', 'traj.phase.parameters:S')\n",
+    "\n",
+    "# Finish Problem Setup\n",
+    "p.setup()\n",
+    "\n",
+    "#############################################\n",
+    "# Set constants and initial guesses\n",
+    "#############################################\n",
+    "p.set_val('radius', 0.05, units='m')\n",
+    "p.set_val('dens', 7.87, units='g/cm**3')\n",
+    "\n",
+    "p.set_val('traj.phase.parameters:CD', 0.5)\n",
+    "\n",
+    "guess = initial_guess(t_dur=30.0)\n",
+    "\n",
+    "p.set_val('traj.phase.t_initial', 0.0)\n",
+    "p.set_val('traj.phase.t_duration', 30.0)\n",
+    "\n",
+    "p.set_val('traj.phase.states:r', phase.interp('r', ys=guess['r'], xs=guess['t']))\n",
+    "p.set_val('traj.phase.states:h', phase.interp('h', ys=guess['h'], xs=guess['t']))\n",
+    "p.set_val('traj.phase.states:v', phase.interp('v', ys=guess['v'], xs=guess['t']))\n",
+    "p.set_val('traj.phase.states:gam', phase.interp('gam', ys=guess['gam'], xs=guess['t']), units='rad')\n",
+    "\n",
+    "#####################################################\n",
+    "# Run the optimization and final explicit simulation\n",
+    "#####################################################\n",
+    "dm.run_problem(p)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Plotting the results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "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.states:gam', units='deg')[0, 0]\n",
+    "print(f'launch angle: {angle} deg')\n",
+    "max_range = p.get_val('traj.phase.timeseries.states: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.states:r')\n",
+    "\n",
+    "r_exp = exp_out.get_val('traj.phase.timeseries.states:r')\n",
+    "\n",
+    "h_imp = p.get_val('traj.phase.timeseries.states:h')\n",
+    "\n",
+    "h_exp = exp_out.get_val('traj.phase.timeseries.states: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.states:{state}')\n",
+    "\n",
+    "    x_exp = exp_out.get_val(f'traj.phase.timeseries.states:{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",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "remove-input",
+     "remove-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from openmdao.utils.assert_utils import assert_near_equal\n",
+    "\n",
+    "assert_near_equal(p.get_val('traj.phase.states:r')[-1],\n",
+    "                  3183.25, tolerance=1.0E-2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "celltoolbar": "Tags",
+  "jupytext": {
+   "cell_metadata_filter": "-all",
+   "notebook_metadata_filter": "-all",
+   "text_representation": {
+    "extension": ".md",
+    "format_name": "markdown"
+   }
+  },
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py
new file mode 100644
index 000000000..f2adecafd
--- /dev/null
+++ b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py
@@ -0,0 +1,294 @@
+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
+    v0 = -0.5*g*t_dur/np.sin(gam0)
+    r = v0*np.cos(gam0)*t
+    h = v0*np.sin(gam0)*t + 0.5*g*t**2
+    v = np.sqrt(v0*np.cos(gam0)**2 + (v0*np.sin(gam0) + g*t)**2)
+    gam = np.arctan2(v0*np.sin(gam0) + g*t, v0*np.cos(gam0)**2)
+
+    guess = {'t': t, 'r': r, 'h': h, 'v': v, 'gam': gam}
+
+    return guess
+
+
+@use_tempdirs
+class TestTwoPhaseCannonballForDocs(unittest.TestCase):
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_two_phase_cannonball_for_docs(self):
+        from scipy.interpolate import interp1d
+
+        import openmdao.api as om
+        from openmdao.utils.assert_utils import assert_near_equal
+
+        import dymos as dm
+        from dymos.models.atmosphere.atmos_1976 import USatm1976Data
+
+        #############################################
+        # Component for the design part of the model
+        #############################################
+        class CannonballSizeComp(om.ExplicitComponent):
+            """
+            Compute the area and mass of a cannonball with a given radius and density.
+
+            Notes
+            -----
+            This component is not vectorized with 'num_nodes' as is the usual way
+            with Dymos, but is instead intended to compute a scalar mass and reference
+            area from scalar radius and density inputs. This component does not reside
+            in the ODE but instead its outputs are connected to the trajectory via
+            input design parameters.
+            """
+            def setup(self):
+                self.add_input(name='radius', val=1.0, desc='cannonball radius', units='m')
+                self.add_input(name='dens', val=7870., desc='cannonball density', units='kg/m**3')
+
+                self.add_output(name='mass', shape=(1,), desc='cannonball mass', units='kg')
+                self.add_output(name='S', shape=(1,), desc='aerodynamic reference area', units='m**2')
+
+                self.declare_partials(of='mass', wrt='dens')
+                self.declare_partials(of='mass', wrt='radius')
+
+                self.declare_partials(of='S', wrt='radius')
+
+            def compute(self, inputs, outputs):
+                radius = inputs['radius']
+                dens = inputs['dens']
+
+                outputs['mass'] = (4/3.) * dens * np.pi * radius ** 3
+                outputs['S'] = np.pi * radius ** 2
+
+            def compute_partials(self, inputs, partials):
+                radius = inputs['radius']
+                dens = inputs['dens']
+
+                partials['mass', 'dens'] = (4/3.) * np.pi * radius ** 3
+                partials['mass', 'radius'] = 4. * dens * np.pi * radius ** 2
+
+                partials['S', 'radius'] = 2 * np.pi * radius
+
+        #############################################
+        # Build the ODE class
+        #############################################
+        class CannonballODE(om.ExplicitComponent):
+            """
+            Cannonball ODE assuming flat earth and accounting for air resistance
+            """
+
+            def initialize(self):
+                self.options.declare('num_nodes', types=int)
+
+            def setup(self):
+                nn = self.options['num_nodes']
+
+                # static parameters
+                self.add_input('m', units='kg')
+                self.add_input('S', units='m**2')
+                # 0.5 good assumption for a sphere
+                self.add_input('CD', 0.5)
+
+                # time varying inputs
+                self.add_input('h', units='m', shape=nn)
+                self.add_input('v', units='m/s', shape=nn)
+                self.add_input('gam', units='rad', shape=nn)
+
+                # state rates
+                self.add_output('v_dot', shape=nn, units='m/s**2', tags=['dymos.state_rate_source:v'])
+                self.add_output('gam_dot', shape=nn, units='rad/s', tags=['dymos.state_rate_source:gam'])
+                self.add_output('h_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:h'])
+                self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r'])
+                self.add_output('ke', shape=nn, units='J')
+
+                # Ask OpenMDAO to compute the partial derivatives using complex-step
+                # with a partial coloring algorithm for improved performance, and use
+                # a graph coloring algorithm to automatically detect the sparsity pattern.
+                self.declare_coloring(wrt='*', method='cs')
+
+                alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]
+                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')
+
+            def compute(self, inputs, outputs):
+
+                gam = inputs['gam']
+                v = inputs['v']
+                h = inputs['h']
+                m = inputs['m']
+                S = inputs['S']
+                CD = inputs['CD']
+
+                GRAVITY = 9.80665  # m/s**2
+
+                # handle complex-step gracefully from the interpolant
+                if np.iscomplexobj(h):
+                    rho = self.rho_interp(inputs['h'])
+                else:
+                    rho = self.rho_interp(inputs['h']).real
+
+                q = 0.5*rho*inputs['v']**2
+                qS = q * S
+                D = qS * CD
+                cgam = np.cos(gam)
+                sgam = np.sin(gam)
+                outputs['v_dot'] = - D/m-GRAVITY*sgam
+                outputs['gam_dot'] = -(GRAVITY/v)*cgam
+                outputs['h_dot'] = v*sgam
+                outputs['r_dot'] = v*cgam
+                outputs['ke'] = 0.5*m*v**2
+
+        #############################################
+        # Setup the Dymos problem
+        #############################################
+
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'IPOPT'
+        p.driver.declare_coloring()
+
+        p.driver.opt_settings['derivative_test'] = 'first-order'
+        p.driver.opt_settings['mu_strategy'] = 'monotone'
+        p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
+        p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
+        p.driver.opt_settings['mu_init'] = 0.01
+        p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
+
+        p.set_solver_print(level=0, depth=99)
+
+        p.model.add_subsystem('size_comp', CannonballSizeComp(),
+                              promotes_inputs=['radius', 'dens'])
+        p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')
+        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())
+
+        transcription = dm.Radau(num_segments=25, order=3, compressed=False)
+        phase = dm.Phase(ode_class=CannonballODE, transcription=transcription)
+
+        phase = traj.add_phase('phase', phase)
+
+        # All initial states except flight path angle are fixed
+        # The output of the ODE which provides the rate source for each state
+        # is obtained from the tags used on those outputs in the ODE.
+        # The units of the states are automatically inferred by multiplying the units
+        # of those rates by the time units.
+        phase.set_time_options(fix_initial=True, duration_bounds=(1, 100), units='s')
+        phase.add_state('r', fix_initial=True, solve_segments='forward')
+        phase.add_state('h', fix_initial=True, solve_segments='forward')
+        phase.add_state('gam', fix_initial=False, solve_segments='forward')
+        phase.add_state('v', fix_initial=False, solve_segments='forward')
+
+        phase.add_parameter('S', units='m**2', static_target=True)
+        phase.add_parameter('m', units='kg', static_target=True)
+        phase.add_parameter('CD', val=0.5, opt=False, static_target=True)
+
+        phase.add_boundary_constraint('ke', loc='initial',
+                                      upper=400000, lower=0, ref=100000)
+
+        # A duration balance is added setting altitude to zero.
+        # A nonlinear solver is used to find the duration of required to satisfy the condition.
+        # The duration was bounded to be greater than 1 to ensure the solver did not
+        # converge to the initial point.
+        phase.set_duration_balance('h', val=0.0)
+
+        phase.add_objective('r', loc='final', scaler=-1.0)
+
+        p.model.connect('size_comp.mass', 'traj.phase.parameters:m')
+        p.model.connect('size_comp.S', 'traj.phase.parameters:S')
+
+        # Finish Problem Setup
+        p.setup()
+
+        #############################################
+        # Set constants and initial guesses
+        #############################################
+        p.set_val('radius', 0.05, units='m')
+        p.set_val('dens', 7.87, units='g/cm**3')
+
+        p.set_val('traj.phase.parameters:CD', 0.5)
+
+        guess = initial_guess(t_dur=30.0)
+
+        p.set_val('traj.phase.t_initial', 0.0)
+        p.set_val('traj.phase.t_duration', 30.0)
+
+        p.set_val('traj.phase.states:r', phase.interp('r', ys=guess['r'], xs=guess['t']))
+        p.set_val('traj.phase.states:h', phase.interp('h', ys=guess['h'], xs=guess['t']))
+        p.set_val('traj.phase.states:v', phase.interp('v', ys=guess['v'], xs=guess['t']))
+        p.set_val('traj.phase.states:gam', phase.interp('gam', ys=guess['gam'], xs=guess['t']), units='rad')
+
+        #####################################################
+        # Run the optimization and final explicit simulation
+        #####################################################
+        dm.run_problem(p)
+
+        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.states:gam', units='deg')[0, 0]
+        print(f'launch angle: {angle} deg')
+        max_range = p.get_val('traj.phase.timeseries.states: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.states:r')
+
+        r_exp = exp_out.get_val('traj.phase.timeseries.states:r')
+
+        h_imp = p.get_val('traj.phase.timeseries.states:h')
+
+        h_exp = exp_out.get_val('traj.phase.timeseries.states: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.states:{state}')
+
+            x_exp = exp_out.get_val(f'traj.phase.timeseries.states:{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/phase/options.py b/dymos/phase/options.py
index 2fc316b87..4d53e0c91 100644
--- a/dymos/phase/options.py
+++ b/dymos/phase/options.py
@@ -571,6 +571,9 @@ def __init__(self, read_only=False):
         self.declare(name='t_duration_targets', allow_none=True, default=[],
                      desc='targets in the ODE to which the total duration of the phase is connected')
 
+        self.declare(name='t_duration_balance_options', default={},
+                     desc='options dictionary for the duration residual')
+
 
 class GridRefinementOptionsDictionary(om.OptionsDictionary):
     """
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 1e9a9c766..2b93461f7 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -710,7 +710,7 @@ def add_polynomial_control(self, name, order, desc=_unspecified, val=_unspecifie
                                ref=_unspecified, targets=_unspecified, rate_targets=_unspecified,
                                rate2_targets=_unspecified, shape=_unspecified):
         """
-        Adds an polynomial control variable to be tied to a parameter in the ODE.
+        Adds a polynomial control variable to be tied to a parameter in the ODE.
 
         Polynomial controls are defined by values at the Legendre-Gauss-Lobatto nodes of a
         single polynomial, defined on [-1, 1] in phase tau space.
@@ -1707,6 +1707,67 @@ def set_time_options(self, units=_unspecified, fix_initial=_unspecified,
         if name is not _unspecified:
             self.time_options['name'] = name
 
+    def set_duration_balance(self, name, val=0.0, index=None, units=None, mult_val=None, normalize=False):
+        """
+        Adds a condition for the duration of the phase. This is satisfied using a nonlinear solver.
+
+        Parameters
+        ----------
+        name : str
+            Name of the variable.  This should be a state or control variable, the path to an output
+            from the top level of the RHS, or an expression to be evaluated. If an expression,
+            it must be provided in the form of an equation with a left- and right-hand side.
+        val :  float
+            The value that the residual must equal at the end  point of the phase.
+        index : int, optional
+            If variable is an array at each point in time, this indicates which index is to be
+            used as the objective, assuming C-ordered flattening.
+        units : str, optional
+            The units of the objective function.  If None, use the units associated with the target.
+            If provided, must be compatible with the target units.
+        mult_val : float, optional
+            Default value for the LHS multiplier.
+        normalize : bool, optional
+            Specifies whether the resulting residual should be normalized by a quadratic
+            function of the RHS.
+        """
+        if self.time_options['fix_duration']:
+            raise ValueError('Cannot implicitly solve for phase duration when fix_duration is True')
+        elif self.time_options['input_duration']:
+            raise ValueError('Cannot implicitly solve for phase duration when input_duration is True')
+
+        if isinstance(self.options['transcription'], ExplicitShooting):
+            raise NotImplementedError('Transcription ExplicitShooting does not implement method setup_duration_balance')
+
+        options = {'name': name,
+                   'val': val,
+                   'index': index,
+                   'units': units,
+                   'mult_val': mult_val,
+                   'normalize': normalize}
+
+        expr_operators = ['(', '+', '-', '/', '*', '&', '%', '@']
+        if '=' in name:
+            is_expr = True
+        elif '=' not in name and any(opr in name for opr in expr_operators):
+            raise ValueError(f'The expression provided `{name}` has invalid format. '
+                             'Expression may be a single variable or an equation '
+                             'of the form `constraint_name = func(vars)`')
+        else:
+            is_expr = False
+
+        balance_name = name.split('=')[0].strip() if is_expr else name
+
+        options['is_expr'] = is_expr
+        options['balance_name'] = balance_name
+
+        self.time_options['t_duration_balance_options'] = options
+
+        var_type = self.classify_var(name)
+        if var_type == 'ode':
+            if balance_name not in self._timeseries['timeseries']['outputs']:
+                self.add_timeseries_output(name, output_name=balance_name, units=units)
+
     def classify_var(self, var):
         """
         Classifies a variable of the given name or path.
@@ -1786,6 +1847,9 @@ def setup(self):
         transcription.setup_ode(self)
 
         transcription.setup_timeseries_outputs(self)
+
+        transcription.setup_duration_balance(self)
+
         transcription.setup_defects(self)
         transcription.setup_solvers(self)
 
@@ -1845,6 +1909,8 @@ def configure(self):
 
         transcription.configure_timeseries_outputs(self)
 
+        transcription.configure_duration_balance(self)
+
         transcription.configure_solvers(self)
 
     def check_time_options(self):
diff --git a/dymos/phase/test/test_implicit_duration.py b/dymos/phase/test/test_implicit_duration.py
new file mode 100644
index 000000000..53c98833f
--- /dev/null
+++ b/dymos/phase/test/test_implicit_duration.py
@@ -0,0 +1,205 @@
+import openmdao.api as om
+import dymos as dm
+import numpy as np
+
+import unittest
+from openmdao.utils.assert_utils import assert_near_equal
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+
+
+class ODEComp(om.ExplicitComponent):
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+        self.add_input('h', shape=nn, units='m')
+        self.add_input('v', shape=nn, units='m/s')
+        self.add_output('hdot', shape=nn, units='m/s')
+        self.add_output('vdot', shape=nn, units='m/s**2')
+
+        self.declare_partials(of='hdot', wrt='v', rows=np.arange(nn), cols=np.arange(nn), val=1.0)
+
+    def compute(self, inputs, outputs):
+        outputs['hdot'] = inputs['v']
+        outputs['vdot'] = -9.80665
+
+
+class MatrixStateCannonball(om.ExplicitComponent):
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+        # z is the state vector, a nn x 2 x 2 in the form of [[x, y], [vx, vy]]
+        self.add_input('z', shape=(nn, 2, 2), units=None)
+        self.add_output('zdot', shape=(nn, 2, 2), units=None)
+
+        self.declare_partials(of='zdot', wrt='z', method='cs')
+        self.declare_coloring(wrt=['z'], method='cs', num_full_jacs=5, tol=1.0E-12)
+
+    def compute(self, inputs, outputs):
+        outputs['zdot'][:, 0, 0] = inputs['z'][:, 1, 0]
+        outputs['zdot'][:, 0, 1] = inputs['z'][:, 1, 1]
+        outputs['zdot'][:, 1, 0] = 0.0
+        outputs['zdot'][:, 1, 1] = -9.81
+
+
+@use_tempdirs
+class TestImplicitDuration(unittest.TestCase):
+
+    @staticmethod
+    def _make_simple_problem(tx, input_duration=False, fix_duration=False, expr_end=False):
+        p = om.Problem()
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=ODEComp, transcription=tx)
+        traj.add_phase('phase', phase)
+        p.model.add_subsystem('traj', traj)
+
+        phase.set_time_options(fix_initial=True, units='s',
+                               input_duration=input_duration, fix_duration=fix_duration)
+        phase.add_state('h', rate_source='hdot', fix_initial=True, units='m', solve_segments='forward')
+        phase.add_state('v', rate_source='vdot', fix_initial=True, units='m/s', solve_segments='forward')
+
+        if not expr_end:
+            phase.set_duration_balance('h', val=0.0, units='m')
+        else:
+            phase.set_duration_balance('pe=9.80665*h', val=0.0, units='m**2/s**2')
+
+        phase.set_simulate_options(rtol=1.0E-9, atol=1.0E-9)
+
+        p.driver = om.ScipyOptimizeDriver(optimizer='SLSQP')
+        p.driver.declare_coloring(tol=1.0E-12)
+
+        p.setup()
+
+        p.set_val('traj.phase.t_initial', 0)
+        p.set_val('traj.phase.t_duration', 2)
+        p.set_val('traj.phase.states:h', phase.interp('h', [30, 0]))
+        p.set_val('traj.phase.states:v', phase.interp('v', [0, -10]))
+
+        return p
+
+    @staticmethod
+    def _make_matrix_problem(tx, shape_error=False):
+        p = om.Problem()
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=MatrixStateCannonball, transcription=tx)
+        traj.add_phase('phase', phase)
+        p.model.add_subsystem('traj', traj)
+
+        phase.set_time_options(fix_initial=True, units='s')
+        phase.add_state('z', rate_source='zdot', fix_initial=True, solve_segments='forward')
+
+        index = None if shape_error else [[0], [1]]
+        phase.set_duration_balance('z', val=0.0, index=index)
+
+        phase.set_simulate_options(rtol=1.0E-9, atol=1.0E-9)
+
+        p.driver = om.ScipyOptimizeDriver(optimizer='SLSQP')
+        p.driver.declare_coloring(tol=1.0E-12)
+
+        p.setup()
+
+        p.set_val('traj.phase.t_initial', 0)
+        p.set_val('traj.phase.t_duration', 7)
+        p.set_val('traj.phase.states:z', phase.interp('z', [[[0, 0], [10, 10]], [[10, 0], [10, -10]]]))
+
+        return p
+
+    def test_implicit_duration_radau(self):
+        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
+
+        p = self._make_simple_problem(tx)
+
+        dm.run_problem(p, run_driver=False, simulate=False)
+
+        assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.4735192, tolerance=1E-6)
+        assert_near_equal(p.get_val('traj.phase.timeseries.states:h')[-1], 0.0, tolerance=1E-6)
+
+    def test_implicit_duration_gl(self):
+        tx = dm.GaussLobatto(num_segments=12, order=3, solve_segments=False)
+
+        p = self._make_simple_problem(tx)
+
+        dm.run_problem(p, run_driver=False, simulate=False)
+
+        assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.4735192, tolerance=1E-6)
+        assert_near_equal(p.get_val('traj.phase.timeseries.states:h')[-1], 0.0, tolerance=1E-6)
+
+    def test_implicit_duration_shooting(self):
+        tx = dm.ExplicitShooting(num_segments=12, order=3)
+
+        expected = 'Transcription ExplicitShooting does not implement method setup_duration_balance'
+
+        with self.assertRaises(NotImplementedError) as e:
+            self._make_simple_problem(tx)
+
+        self.assertEqual(expected, str(e.exception))
+
+    def test_implicit_duration_input_duration(self):
+        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
+
+        expected = 'Cannot implicitly solve for phase duration when input_duration is True'
+
+        with self.assertRaises(ValueError) as e:
+            self._make_simple_problem(tx, input_duration=True)
+
+        self.assertEqual(expected, str(e.exception))
+
+    def test_implicit_duration_fix_duration(self):
+        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
+
+        expected = 'Cannot implicitly solve for phase duration when fix_duration is True'
+
+        with self.assertRaises(ValueError) as e:
+            self._make_simple_problem(tx, fix_duration=True)
+
+        self.assertEqual(expected, str(e.exception))
+
+    def test_implicit_duration_matrix_state_exception(self):
+        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
+
+        expected = 'Specified variable for the duration balance has shape (2, 2) and' \
+                   ' has no index specified. The balance may only have shape (1,) or a single index'
+
+        with self.assertRaises(ValueError) as e:
+            self._make_matrix_problem(tx, shape_error=True)
+
+        self.assertEqual(expected, str(e.exception))
+
+    def test_implicit_duration_matrix_state(self):
+        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
+
+        p = self._make_matrix_problem(tx, shape_error=False)
+
+        dm.run_problem(p, run_driver=False, simulate=False)
+
+        assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_implicit_duration_radau_expr_condition(self):
+        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
+
+        p = self._make_simple_problem(tx, expr_end=True)
+
+        dm.run_problem(p, run_driver=False, simulate=False)
+
+        assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.4735192, tolerance=1E-6)
+        assert_near_equal(p.get_val('traj.phase.timeseries.states:h')[-1], 0.0, tolerance=1E-6)
+
+        print((p.get_val('traj.phase.timeseries.states:v')[-1]) ** 2 / 2)
+
+    def test_implicit_duration_gl_expr_condition(self):
+        tx = dm.GaussLobatto(num_segments=12, order=3, solve_segments=False)
+
+        p = self._make_simple_problem(tx, expr_end=True)
+
+        dm.run_problem(p, run_driver=False, simulate=False)
+
+        assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.4735192, tolerance=1E-6)
+        assert_near_equal(p.get_val('traj.phase.timeseries.states:h')[-1], 0.0, tolerance=1E-6)
diff --git a/dymos/transcriptions/analytic/analytic.py b/dymos/transcriptions/analytic/analytic.py
index 8b5f55240..76493a794 100644
--- a/dymos/transcriptions/analytic/analytic.py
+++ b/dymos/transcriptions/analytic/analytic.py
@@ -247,6 +247,28 @@ def setup_defects(self, phase):
         """
         pass
 
+    def setup_duration_balance(self, phase):
+        """
+        Setup the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
+    def configure_duration_balance(self, phase):
+        """
+        Configure the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
     def setup_solvers(self, phase):
         """
         Setup the solvers for this transcription.
@@ -260,7 +282,7 @@ def setup_solvers(self, phase):
 
     def configure_solvers(self, phase):
         """
-        Setup the solvers for this transcription.
+        Configure the solvers for this transcription.
 
         Parameters
         ----------
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index 8fdba24c5..a22fcab3a 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -643,6 +643,28 @@ def configure_timeseries_outputs(self, phase):
         """
         super().configure_timeseries_outputs(phase)
 
+    def setup_duration_balance(self, phase):
+        """
+        Setup the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
+    def configure_duration_balance(self, phase):
+        """
+        Configure the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
     def setup_solvers(self, phase):
         """
         Setup the solvers for this transcription.
diff --git a/dymos/transcriptions/pseudospectral/components/state_independents.py b/dymos/transcriptions/pseudospectral/components/state_independents.py
index 45a8ffb3f..2c3e34602 100644
--- a/dymos/transcriptions/pseudospectral/components/state_independents.py
+++ b/dymos/transcriptions/pseudospectral/components/state_independents.py
@@ -12,7 +12,7 @@ class StateIndependentsComp(om.ImplicitComponent):
     """
     Class definition for the StateIndependentsComp.
 
-    A simple component that replaces the state indepvarcomps whenver the solver needs to solve for
+    A simple component that replaces the state IndepVarComps whenever the solver needs to solve for
     the state or whenever the initial state is connected to an external source.
 
     Parameters
diff --git a/dymos/transcriptions/pseudospectral/components/time_duration_comp.py b/dymos/transcriptions/pseudospectral/components/time_duration_comp.py
new file mode 100644
index 000000000..50041c586
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/time_duration_comp.py
@@ -0,0 +1,211 @@
+import numpy as np
+import openmdao.api as om
+
+from dymos.transcriptions.grid_data import GridData
+from dymos._options import options as dymos_options
+
+
+class TimeDurationComp(om.ImplicitComponent):
+    """
+    Class definition for the TimeDurationComp.
+
+    A simple component that computes the duration of the phase when the solver needs to solve for
+    the stopping condition.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._no_check_partials = not dymos_options['include_check_partials']
+
+    def initialize(self):
+        """
+        Declare component options.
+        """
+        self.options.declare(
+            'grid_data', types=GridData,
+            desc='Container object for grid info')
+
+        self.options.declare(
+            'state_options', types=dict,
+            desc='Dictionary of state names/options for the phase')
+
+    def configure_io(self, state_idx_map=None):
+        """
+        I/O creation is delayed until configure so we can determine shape and units.
+
+        Parameters
+        ----------
+        state_idx_map : dict
+            The state_idx_map generated by the transcription during its configuration.  This dict
+            is keyed by state, which in turn holds a dict of the solve indices in the input state
+            vector (keyed as 'solved') and those indices managed by the optimizer (keyed as 'indep').
+            In some unittest, this map is not provided at the time of setup and is instead called
+            later.
+        """
+        if state_idx_map is None:
+            return
+        self.state_idx_map = state_idx_map
+        state_options = self.options['state_options']
+        grid_data = self.options['grid_data']
+        num_col_nodes = grid_data.subset_num_nodes['col']
+
+        num_state_input_nodes = grid_data.subset_num_nodes['state_input']
+
+        self.var_names = {}
+        for state_name in state_options:
+            self.var_names[state_name] = {
+                'defect': f'defects:{state_name}',
+            }
+
+        for state_name, options in state_options.items():
+
+            shape = options['shape']
+            units = options['units']
+            solved = options['solve_segments']
+            default_val = options['val']
+            if np.isscalar(default_val) or np.asarray(default_val).shape == 1:
+                default_val = float(default_val)
+            elif np.asarray(default_val).shape == shape:
+                default_val = np.repeat(default_val[np.newaxis, ...], num_state_input_nodes, axis=0)
+            var_names = self.var_names[state_name]
+
+            if solved and (options['lower'] or options['upper']):
+                om.issue_warning(f'State {state_name} has bounds but they are not enforced when '
+                                 f'using `solve_segments.` Apply a path constraint to {state_name} '
+                                 f'to enforce bounds.', om.UnusedOptionWarning)
+
+            # only need the implicit variable if this state is solved.
+            # Note: we don't add scaling and bounds here. This may be revisited.
+            self.add_output(name=f'states:{state_name}',
+                            shape=(num_state_input_nodes,) + shape,
+                            val=default_val,
+                            units=units)
+
+            # Input for continuity, which can come from an external source.
+            if options['input_initial']:
+                input_name = f'initial_states:{state_name}'
+                self.add_input(name=input_name, shape=(1, ) + shape, units=units)
+
+            # compute an output constraint value since the optimizer needs it
+            if solved:
+                self.add_input(
+                    name=var_names['defect'],
+                    shape=(num_col_nodes, ) + shape,
+                    desc=f'Constraint value for interior defects of state {state_name}',
+                    units=units)
+
+        # Setup partials
+        for state_name, options in state_options.items():
+            shape = options['shape']
+            size = np.prod(shape)
+            solved = options['solve_segments']
+            state_var_name = f'states:{state_name}'
+
+            if solved:  # only need this deriv if its solved
+                solve_idx = np.array(state_idx_map[state_name]['solver'])
+                indep_idx = np.array(state_idx_map[state_name]['indep'])
+
+                num_indep_nodes = indep_idx.shape[0]
+                num_solve_nodes = solve_idx.shape[0]
+
+                base_idx = np.tile(np.arange(size), num_indep_nodes).reshape(num_indep_nodes, size)
+                row = (indep_idx[:, np.newaxis]*size + base_idx).flatten()
+
+                # anything that looks like an indep
+                self.declare_partials(of=state_var_name, wrt=state_var_name,
+                                      rows=row, cols=row, val=-1.0)
+
+                if options['input_initial']:
+                    wrt = f'initial_states:{state_name}'
+                    row_col = np.arange(np.prod(shape))
+                    self.declare_partials(of=state_var_name, wrt=wrt, rows=row_col, cols=row_col,
+                                          val=1.0)
+
+                col = np.arange(num_solve_nodes * size)
+                base_idx = np.tile(np.arange(size), num_solve_nodes).reshape(num_solve_nodes, size)
+                row = (solve_idx[:, np.newaxis]*size + base_idx).flatten()
+
+                var_names = self.var_names[state_name]
+                self.declare_partials(of=state_var_name,
+                                      wrt=var_names['defect'],
+                                      rows=row, cols=col, val=1.0)
+
+            else:
+                row_col = np.arange(num_state_input_nodes*np.prod(shape), dtype=int)
+                self.declare_partials(of=state_var_name, wrt=state_var_name,
+                                      rows=row_col, cols=row_col, val=-1.0)
+
+                if options['input_initial']:
+                    wrt = f'initial_states:{state_name}'
+                    row_col = np.arange(np.prod(shape))
+                    self.declare_partials(of=state_var_name, wrt=wrt, rows=row_col, cols=row_col,
+                                          val=1.0)
+
+    def apply_nonlinear(self, inputs, outputs, residuals):
+        """
+        Compute residuals given inputs and outputs.
+
+        The model is assumed to be in an unscaled state.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        outputs : Vector
+            Unscaled, dimensional output variables read via outputs[key].
+        residuals : Vector
+            Unscaled, dimensional residuals written to via residuals[key].
+        """
+        state_options = self.options['state_options']
+
+        for state_name, options in state_options.items():
+            state_var_name = f'states:{state_name}'
+
+            solve_idx = self.state_idx_map[state_name]['solver']
+            indep_idx = self.state_idx_map[state_name]['indep']
+
+            var_names = self.var_names[state_name]
+
+            if options['solve_segments']:
+                defect = inputs[var_names['defect']]
+                residuals[state_var_name][solve_idx, ...] = defect
+
+            # really is: <idep_val> - \outputs[state_name][indep_idx] but OpenMDAO
+            # implementation details mean we just set it to 0
+            # but derivatives are still based on (<idep_val> - \outputs[state_name][indep_idx]),
+            # so you get -1 wrt state var
+            # NOTE: check_partials will report wrong derivs for the indep vars,
+            #       but don't believe it!
+            residuals[state_var_name][indep_idx, ...] = 0.0
+
+            if options['input_initial']:
+                ic_state_name = f'initial_states:{state_name}'
+
+                residuals[state_var_name][0, ...] = \
+                    inputs[ic_state_name][0, ...] - \
+                    outputs[state_var_name][0, ...]
+
+    def solve_nonlinear(self, inputs, outputs):
+        """
+        Compute outputs given inputs.
+
+        The model is assumed to be in an unscaled state.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        outputs : Vector
+            Unscaled, dimensional output variables read via outputs[key].
+        """
+        state_options = self.options['state_options']
+
+        for state_name, options in state_options.items():
+            if options['input_initial']:
+                output_name = f'states:{state_name}'
+                input_name = f'initial_states:{state_name}'
+                outputs[output_name][0, ...] = inputs[input_name]
diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
index ba3db3bff..86ed875f7 100644
--- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py
+++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
@@ -6,7 +6,7 @@
 from .components import StateIndependentsComp, StateInterpComp, CollocationComp, \
     PseudospectralTimeseriesOutputComp
 from ...utils.misc import CoerceDesvar, get_rate_units, reshape_val
-from ...utils.introspection import get_promoted_vars, get_source_metadata
+from ...utils.introspection import get_promoted_vars, get_source_metadata, configure_duration_balance_introspection
 from ...utils.constants import INF_BOUND
 from ...utils.indexing import get_src_indices_by_row
 
@@ -436,6 +436,53 @@ def configure_defects(self, phase):
                               f'continuity_comp.control_rates:{name}_rate2',
                               src_indices=src_idxs, flat_src_indices=True)
 
+    def setup_duration_balance(self, phase):
+        """
+        Setup the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+
+        if self._implicit_duration:
+            duration_balance_comp = om.BalanceComp()
+            phase.add_subsystem('t_duration_balance_comp', duration_balance_comp)
+
+    def configure_duration_balance(self, phase):
+        """
+        Configure the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+
+        if self._implicit_duration:
+            duration_balance_comp = phase._get_subsystem('t_duration_balance_comp')
+            configure_duration_balance_introspection(phase)
+            options = phase.time_options['t_duration_balance_options']
+
+            if options['mult_val'] is None:
+                use_mult = False
+                mult_val = 1.0
+            else:
+                use_mult = True
+                mult_val = options['mult_val']
+
+            src_idx = [-1] if options['index'] is None else [[-1]]+options['index']
+
+            duration_balance_comp.add_balance('t_duration', val=options['val'],
+                                              eq_units=options['units'], units=phase.time_options['units'],
+                                              use_mult=use_mult, mult_val=mult_val)
+
+            phase.connect('t_duration_balance_comp.t_duration', 't_duration')
+
+            phase.connect(options['var_path'], 't_duration_balance_comp.lhs:t_duration',
+                          src_indices=src_idx)
+
     def setup_solvers(self, phase):
         """
         Setup the solvers.
@@ -456,7 +503,7 @@ def configure_solvers(self, phase):
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
         """
-        if self.any_solved_segs:
+        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()
@@ -467,7 +514,7 @@ def configure_solvers(self, phase):
                 newton.linesearch = om.BoundsEnforceLS()
 
         # even though you don't need a nl_solver for connections, you still ln_solver since its implicit
-        if self.any_solved_segs or self.any_connected_opt_segs:
+        if self.any_solved_segs or self.any_connected_opt_segs or self._implicit_duration:
             if isinstance(phase.linear_solver, om.LinearRunOnce):
                 phase.linear_solver = om.DirectSolver()
 
diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index 360b090d0..68123c66e 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -24,6 +24,7 @@ class TranscriptionBase(object):
     """
     def __init__(self, **kwargs):
 
+        self._implicit_duration = False
         self.grid_data = None
 
         self.options = om.OptionsDictionary()
@@ -80,9 +81,16 @@ def setup_time(self, phase):
         # Warn about invalid options
         phase.check_time_options()
 
-        if not time_options['input_initial'] or not time_options['input_duration']:
+        if phase.time_options['t_duration_balance_options']:
+            self._implicit_duration = True
+
+        if not time_options['input_initial']:
             phase.add_subsystem('time_extents', om.IndepVarComp(),
                                 promotes_outputs=['*'])
+        else:
+            if not time_options['input_duration'] and not self._implicit_duration:
+                phase.add_subsystem('time_extents', om.IndepVarComp(),
+                                    promotes_outputs=['*'])
 
         for ts_name, ts_options in phase._timeseries.items():
             if t_name not in ts_options['outputs']:
@@ -120,7 +128,7 @@ def configure_time(self, phase):
         if not time_options['input_initial']:
             indeps.append('t_initial')
 
-        if not time_options['input_duration']:
+        if not time_options['input_duration'] and not self._implicit_duration:
             indeps.append('t_duration')
 
         for var in indeps:
@@ -379,6 +387,31 @@ def setup_ode(self, phase):
         """
         raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement method setup_ode.')
 
+    def setup_duration_balance(self, phase):
+        """
+        Setup the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+
+        raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement'
+                                  f' method setup_duration_balance.')
+
+    def configure_duration_balance(self, phase):
+        """
+        Configure the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement'
+                                  f' method setup_duration_balance.')
+
     def setup_solvers(self, phase):
         """
         Setup the solvers for this transcription.
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index 1bc1776af..0120e660b 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -1061,6 +1061,122 @@ def get_target_metadata(ode, name, user_targets=_unspecified, user_units=_unspec
     return shape, units, static_target
 
 
+def configure_duration_balance_introspection(phase):
+    """
+    Modify duration balance options in-place using introspection of the phase and its ODE.
+
+    Parameters
+    ----------
+    phase : Phase
+        The phase object whose boundary and path constraints are to be introspected.
+    """
+    time_units = phase.time_options['units']
+
+    options = phase.time_options['t_duration_balance_options']
+
+    # Determine the path to the variable which we will be constraining
+    var = options['balance_name'] if options['is_expr'] else options['name']
+    var_type = phase.classify_var(var)
+
+    if var != options['balance_name'] is not None and var_type != 'ode':
+        om.issue_warning(f"Option 'balance_name' on duration residual {var} is only "
+                         f"valid for ODE outputs. The option is being ignored.", om.UnusedOptionWarning)
+
+    if var_type == 't' or var_type == 't_phase':
+        raise ValueError(f'Cannot use a variable of type {var_type} in the duration balance')
+
+    elif var_type == 'state':
+        prefix = 'states:' if dymos_options['use_timeseries_prefix'] else ''
+        state_shape = phase.state_options[var]['shape']
+        state_units = phase.state_options[var]['units']
+        options['shape'] = state_shape
+        options['units'] = state_units if options['units'] is None else options['units']
+        options['var_path'] = f'timeseries.{prefix}{var}'
+
+    elif var_type == 'parameter':
+        param_shape = phase.parameter_options[var]['shape']
+        param_units = phase.parameter_options[var]['units']
+        options['shape'] = param_shape
+        options['units'] = param_units if options['units'] is None else options['units']
+        options['var_path'] = f'parameter_vals:{var}'
+
+    elif var_type in ['indep_control', 'input_control']:
+        prefix = 'controls:' if dymos_options['use_timeseries_prefix'] else ''
+        control_shape = phase.control_options[var]['shape']
+        control_units = phase.control_options[var]['units']
+
+        options['shape'] = control_shape
+        options['units'] = control_units if options['units'] is None else options['units']
+        options['var_path'] = f'timeseries.{prefix}{var}'
+
+    elif var_type in ['indep_polynomial_control', 'input_polynomial_control']:
+        prefix = 'polynomial_controls:' if dymos_options['use_timeseries_prefix'] else ''
+        control_shape = phase.polynomial_control_options[var]['shape']
+        control_units = phase.polynomial_control_options[var]['units']
+        options['shape'] = control_shape
+        options['units'] = control_units if options['units'] is None else options['units']
+        options['var_path'] = f'timeseries.{prefix}{var}'
+
+    elif var_type == 'control_rate':
+        prefix = 'control_rates:' if dymos_options['use_timeseries_prefix'] else ''
+        control_name = var[:-5]
+        control_shape = phase.control_options[control_name]['shape']
+        control_units = phase.control_options[control_name]['units']
+        options['shape'] = control_shape
+        options['units'] = get_rate_units(control_units, time_units, deriv=1) \
+            if options['units'] is None else options['units']
+        options['var_path'] = f'timeseries.{prefix}{var}'
+
+    elif var_type == 'control_rate2':
+        prefix = 'control_rates:' if dymos_options['use_timeseries_prefix'] else ''
+        control_name = var[:-6]
+        control_shape = phase.control_options[control_name]['shape']
+        control_units = phase.control_options[control_name]['units']
+        options['shape'] = control_shape
+        options['units'] = get_rate_units(control_units, time_units, deriv=2) \
+            if options['units'] is None else options['units']
+        options['var_path'] = f'timeseries.{prefix}{var}'
+
+    elif var_type == 'polynomial_control_rate':
+        prefix = 'polynomial_control_rates:' if dymos_options['use_timeseries_prefix'] else ''
+        control_name = var[:-5]
+        control_shape = phase.polynomial_control_options[control_name]['shape']
+        control_units = phase.polynomial_control_options[control_name]['units']
+        options['shape'] = control_shape
+        options['units'] = get_rate_units(control_units, time_units, deriv=1) \
+            if options['units'] is None else options['units']
+        options['var_path'] = f'timeseries.{prefix}{var}'
+
+    elif var_type == 'polynomial_control_rate2':
+        prefix = 'polynomial_control_rates:' if dymos_options['use_timeseries_prefix'] else ''
+        control_name = var[:-6]
+        control_shape = phase.polynomial_control_options[control_name]['shape']
+        control_units = phase.polynomial_control_options[control_name]['units']
+        options['shape'] = control_shape
+        options['units'] = get_rate_units(control_units, time_units, deriv=2) \
+            if options['units'] is None else options['units']
+        options['var_path'] = f'timeseries.{prefix}{var}'
+
+    elif var_type == 'timeseries_exec_comp_output':
+        options['shape'] = (1,)
+        options['units'] = None
+        options['var_path'] = f'timeseries.timeseries_exec_comp.{var}'
+
+    else:
+        # Failed to find variable, assume it is in the ODE. This requires introspection.
+        ode = phase.options['transcription']._get_ode(phase)
+
+        meta = get_source_metadata(ode, src=var, user_units=options['units'])
+
+        options['shape'] = meta['shape']
+        options['units'] = meta['units']
+        options['var_path'] = f'timeseries.{options["balance_name"]}'
+
+    if options['shape'] != (1,) and options['index'] is None:
+        raise ValueError(f'Specified variable for the duration balance has shape {options["shape"]} and'
+                         f' has no index specified. The balance may only have shape (1,) or a single index')
+
+
 def _get_targets_metadata(ode, name, user_targets=_unspecified, user_units=_unspecified,
                           user_shape=_unspecified, control_rate=False, user_static_target=_unspecified):
     """