diff --git a/docs/source/examples/shock/example7.rst b/docs/source/examples/shock/example7.rst
index b5936ff..3a24d7c 100644
--- a/docs/source/examples/shock/example7.rst
+++ b/docs/source/examples/shock/example7.rst
@@ -112,6 +112,9 @@ Initialize arrays for storing the solution variables:
     mole_fractions_refl = {}
 
 Loop over the initial shock velocities and solve the shock problem, noting again that we set `reflected=True` to compute the equilibrium reflected shock solution.
+This example still gates downstream use on ``solution.converged``. If a solve
+returns ``solution.last_error == cea.LAST_VALID_SOLUTION``, the retained state
+is still not a converged shock point.
 
 .. code-block:: python
     :emphasize-lines: 4
diff --git a/docs/source/interfaces/c_api.rst b/docs/source/interfaces/c_api.rst
index ba2d777..1912ad5 100644
--- a/docs/source/interfaces/c_api.rst
+++ b/docs/source/interfaces/c_api.rst
@@ -11,6 +11,9 @@ Notes:
 - ``cea_rocket_solver_get_eqsolver`` is no longer exposed; use ``cea_rocket_solver_get_size`` and other RocketSolver APIs.
 - Rocket, shock, and detonation solution property getters now require a ``len`` argument and will return
   ``CEA_INVALID_SIZE`` if ``len`` is smaller than the internal number of points.
+- ``cea_shock_solver_solve`` may return ``CEA_LAST_VALID_SOLUTION`` when the incident-equilibrium solver retains
+  a last valid state for inspection without converging the shock iteration. In that case the returned properties
+  remain readable, but ``cea_shock_solution_get_converged`` still reports false.
 
 .. doxygenfile:: cea.h
    :project: cea
diff --git a/docs/source/interfaces/python_api.rst b/docs/source/interfaces/python_api.rst
index 84b723c..4946d4b 100644
--- a/docs/source/interfaces/python_api.rst
+++ b/docs/source/interfaces/python_api.rst
@@ -49,6 +49,12 @@ RocketSolution
 ShockSolver
 -----------
 
+Shock solves can emit a ``RuntimeWarning`` with ``LAST_VALID_SOLUTION`` when the
+incident-equilibrium path retains a last valid state without converging the
+shock iteration. In that case ``ShockSolution.last_error`` records the soft
+failure, while ``ShockSolution.converged`` remains the authoritative strict
+convergence signal.
+
 .. autoclass:: cea.ShockSolver
    :members:
 
diff --git a/source/bind/c/CMakeLists.txt b/source/bind/c/CMakeLists.txt
index 5a0a981..d49003e 100644
--- a/source/bind/c/CMakeLists.txt
+++ b/source/bind/c/CMakeLists.txt
@@ -97,4 +97,12 @@ if (CEA_BUILD_TESTING)
         WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
     )
 
+    add_executable(cea_bindc_shock_last_valid samples/shock_last_valid.c)
+    target_link_libraries(cea_bindc_shock_last_valid PRIVATE cea::bindc)
+    add_test(
+        NAME cea_bindc_shock_last_valid
+        COMMAND cea_bindc_shock_last_valid
+        WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+    )
+
 endif()
diff --git a/source/bind/c/bindc.F90 b/source/bind/c/bindc.F90
index bcd62a0..6f9c0ab 100644
--- a/source/bind/c/bindc.F90
+++ b/source/bind/c/bindc.F90
@@ -7,6 +7,7 @@ module cea_bindc
     use cea_param, only: empty_dp, gas_constant, get_data_search_dirs
     use cea_input, only: ReactantInput
     use cea_mixture, only: names_match
+    use cea_shock, only: SHOCK_SOLVE_STATUS_LAST_VALID
     use iso_c_binding
     use fb_logging
     use fb_utils, only: assert, locate, is_empty, to_str
@@ -203,6 +204,7 @@ module cea_bindc
         enumerator :: CEA_INVALID_INDEX = 6
         enumerator :: CEA_INVALID_SIZE = 7
         enumerator :: CEA_NOT_CONVERGED = 8
+        enumerator :: CEA_LAST_VALID_SOLUTION = 9
     end enum
 
     !-----------------------------------------------------------------
@@ -2495,7 +2497,16 @@ function cea_shock_solver_solve(sptr, slptr, weights, T0, p0, mach1_or_u1, use_m
             solution = solver%solve(weights(:nr), T0, p0, u1=mach1_or_u1, &
                 reflected=reflected, incident_frozen=incident_frozen, reflected_frozen=reflected_frozen)
         end if
-        if (.not. solution%converged) ierr = CEA_NOT_CONVERGED
+        if (.not. solution%converged) then
+            ! A retained last-valid incident state is intentionally exposed as a
+            ! descriptive soft failure, but it remains a non-converged shock
+            ! solution and callers must check the converged flag separately.
+            if (solution%solve_status == SHOCK_SOLVE_STATUS_LAST_VALID) then
+                ierr = CEA_LAST_VALID_SOLUTION
+            else
+                ierr = CEA_NOT_CONVERGED
+            end if
+        end if
     end function
 
 
diff --git a/source/bind/c/cea.h b/source/bind/c/cea.h
index 0be18ae..30fb558 100644
--- a/source/bind/c/cea.h
+++ b/source/bind/c/cea.h
@@ -639,6 +639,9 @@ extern "C"
       cea_int *value);
 
   // Solve
+  // Returns CEA_LAST_VALID_SOLUTION when the incident-equilibrium shock path
+  // retains a last valid state for inspection without converging the shock
+  // iteration. In that case, cea_shock_solution_get_converged still reports 0.
   cea_err cea_shock_solver_solve(
       const cea_shock_solver solver,
       cea_shock_solution soln,
diff --git a/source/bind/c/cea_enum.h b/source/bind/c/cea_enum.h
index 7f3a0e4..94cd8f7 100644
--- a/source/bind/c/cea_enum.h
+++ b/source/bind/c/cea_enum.h
@@ -189,4 +189,5 @@
     CEA_INVALID_EQUILIBRIUM_SIZE_TYPE = 5, \
     CEA_INVALID_INDEX            = 6, \
     CEA_INVALID_SIZE             = 7, \
-    CEA_NOT_CONVERGED            = 8
+    CEA_NOT_CONVERGED            = 8, \
+    CEA_LAST_VALID_SOLUTION      = 9
diff --git a/source/bind/c/samples/shock_last_valid.c b/source/bind/c/samples/shock_last_valid.c
new file mode 100644
index 0000000..f5528d2
--- /dev/null
+++ b/source/bind/c/samples/shock_last_valid.c
@@ -0,0 +1,143 @@
+#include "stdlib.h"
+#include "stdio.h"
+#include "cea.h"
+
+#define LEN(x) (sizeof(x) / sizeof((x)[0]))
+#define NUM_PTS 2
+#define NUM_CANDIDATES 12
+
+int main(void) {
+    const cea_string reactants[] = {"N2 ", "CH4"};
+    const cea_real moles[] = {0.985, 0.015};
+    const cea_string omitted_products[] = {""};
+    const cea_int nr = LEN(reactants);
+    const cea_real p0 = 0.0002986421052631579;
+    const cea_real T0 = 291.0;
+    const cea_real u1_candidates[NUM_CANDIDATES] = {
+        1200.0, 1300.0, 1400.0, 1500.0, 1600.0, 1700.0,
+        1800.0, 1900.0, 2000.0, 2100.0, 2200.0, 2300.0
+    };
+
+    int status = EXIT_FAILURE;
+    int found_last_valid = 0;
+    cea_err ierr;
+    cea_mixture reac = NULL;
+    cea_mixture prod = NULL;
+    cea_shock_solver solver = NULL;
+    cea_shock_solution soln = NULL;
+    cea_real *weights = NULL;
+    cea_real pressure[NUM_PTS];
+    cea_real temperature[NUM_PTS];
+    int converged = 1;
+    size_t i;
+    cea_solver_opts opts;
+
+    cea_set_log_level(CEA_LOG_NONE);
+    ierr = cea_init();
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "cea_init failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    ierr = cea_mixture_create_w_ions(&reac, nr, reactants);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "reactant create failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    ierr = cea_mixture_create_from_reactants_w_ions(&prod, nr, reactants, 0, omitted_products);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "product create failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    ierr = cea_solver_opts_init(&opts);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "opts init failed: %d\n", ierr);
+        goto cleanup;
+    }
+    opts.trace = 1.0e-10;
+    opts.ions = true;
+    opts.transport = true;
+    opts.reactants = reac;
+
+    ierr = cea_shock_solver_create_with_options(&solver, prod, opts);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "solver create failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    ierr = cea_shock_solution_create(&soln, NUM_PTS);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "solution create failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    weights = (cea_real *) calloc(nr, sizeof(cea_real));
+    if (weights == NULL) {
+        fprintf(stderr, "weights allocation failed\n");
+        goto cleanup;
+    }
+
+    ierr = cea_mixture_moles_to_weights(reac, nr, moles, weights);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "weights conversion failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    for (i = 0; i < NUM_CANDIDATES; ++i) {
+        ierr = cea_shock_solver_solve(
+            solver, soln, weights, T0, p0, u1_candidates[i], false, false, false, false
+        );
+        if (ierr == CEA_LAST_VALID_SOLUTION) {
+            found_last_valid = 1;
+            break;
+        }
+        if (ierr != CEA_SUCCESS && ierr != CEA_NOT_CONVERGED) {
+            fprintf(stderr, "unexpected solve status for u1=%g: %d\n", u1_candidates[i], ierr);
+            goto cleanup;
+        }
+    }
+
+    if (!found_last_valid) {
+        fprintf(stderr, "no last-valid incident shock case found across candidate velocities\n");
+        goto cleanup;
+    }
+
+    ierr = cea_shock_solution_get_converged(soln, &converged);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "get_converged failed: %d\n", ierr);
+        goto cleanup;
+    }
+    if (converged != 0) {
+        fprintf(stderr, "expected non-converged shock status\n");
+        goto cleanup;
+    }
+
+    ierr = cea_shock_solution_get_property(soln, CEA_SHOCK_PRESSURE, NUM_PTS, pressure);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "pressure query failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    ierr = cea_shock_solution_get_property(soln, CEA_SHOCK_TEMPERATURE, NUM_PTS, temperature);
+    if (ierr != CEA_SUCCESS) {
+        fprintf(stderr, "temperature query failed: %d\n", ierr);
+        goto cleanup;
+    }
+
+    if (pressure[1] <= 0.0 || temperature[1] <= 0.0) {
+        fprintf(stderr, "expected retained incident state to remain populated\n");
+        goto cleanup;
+    }
+
+    status = EXIT_SUCCESS;
+
+cleanup:
+    if (soln != NULL) cea_shock_solution_destroy(&soln);
+    if (solver != NULL) cea_shock_solver_destroy(&solver);
+    if (prod != NULL) cea_mixture_destroy(&prod);
+    if (reac != NULL) cea_mixture_destroy(&reac);
+    free(weights);
+    return status;
+}
diff --git a/source/bind/python/CEA.pyx b/source/bind/python/CEA.pyx
index 0089fab..011922b 100644
--- a/source/bind/python/CEA.pyx
+++ b/source/bind/python/CEA.pyx
@@ -62,10 +62,19 @@ def _err_name(cea_err ierr):
         return "CEA_INVALID_SIZE"
     if ierr == CEA_NOT_CONVERGED:
         return "CEA_NOT_CONVERGED"
+    if ierr == CEA_LAST_VALID_SOLUTION:
+        return "CEA_LAST_VALID_SOLUTION"
     return f"CEA_ERR_{int(ierr)}"
 
 
 def _check_ierr(cea_err ierr, str context, bint allow_not_converged=True):
+    if ierr == CEA_LAST_VALID_SOLUTION and allow_not_converged:
+        warnings.warn(
+            f"{context}: CEA_LAST_VALID_SOLUTION "
+            f"(last valid shock state retained; solution.converged is False)",
+            RuntimeWarning,
+        )
+        return
     if ierr == CEA_NOT_CONVERGED and allow_not_converged:
         warnings.warn(f"{context}: { _err_name(ierr) }", RuntimeWarning)
         return
@@ -77,6 +86,8 @@ SUCCESS                  = CEA_SUCCESS
 INVALID_FILENAME         = CEA_INVALID_FILENAME
 INVALID_PROPERTY_TYPE    = CEA_INVALID_PROPERTY_TYPE
 INVALID_EQUILIBRIUM_TYPE = CEA_INVALID_EQUILIBRIUM_TYPE
+NOT_CONVERGED           = CEA_NOT_CONVERGED
+LAST_VALID_SOLUTION     = CEA_LAST_VALID_SOLUTION
 
 # Alias the log levels
 LOG_CRITICAL = CEA_LOG_CRITICAL
@@ -3335,6 +3346,12 @@ cdef class ShockSolver:
         reflected_frozen : bool, default False
             Use frozen chemistry for reflected shock
 
+        Notes
+        -----
+        A retained incident last-valid state emits ``RuntimeWarning`` and sets
+        ``soln.last_error`` to ``LAST_VALID_SOLUTION`` while leaving
+        ``soln.converged`` as ``False``.
+
         Raises
         ------
         ValueError
@@ -3921,7 +3938,7 @@ cdef class ShockSolution:
         Returns
         -------
         bool
-            True if solution converged, False otherwise
+            True if the shock iteration converged, False otherwise
         """
         def __get__(self):
             cdef cea_err ierr
diff --git a/source/bind/python/cea/lib/libcea.pyi b/source/bind/python/cea/lib/libcea.pyi
index a69b03a..6806fc5 100644
--- a/source/bind/python/cea/lib/libcea.pyi
+++ b/source/bind/python/cea/lib/libcea.pyi
@@ -36,6 +36,8 @@ SUCCESS: int
 INVALID_FILENAME: int
 INVALID_PROPERTY_TYPE: int
 INVALID_EQUILIBRIUM_TYPE: int
+NOT_CONVERGED: int
+LAST_VALID_SOLUTION: int
 
 LOG_CRITICAL: int
 LOG_ERROR: int
diff --git a/source/bind/python/cea_def.pxd b/source/bind/python/cea_def.pxd
index 3d24cc8..bdb756f 100644
--- a/source/bind/python/cea_def.pxd
+++ b/source/bind/python/cea_def.pxd
@@ -13,6 +13,7 @@ cdef extern from "cea.h":
         CEA_INVALID_INDEX
         CEA_INVALID_SIZE
         CEA_NOT_CONVERGED
+        CEA_LAST_VALID_SOLUTION
 
     ctypedef enum cea_log_level:
         CEA_LOG_CRITICAL
diff --git a/source/bind/python/tests/test_shock_last_valid.py b/source/bind/python/tests/test_shock_last_valid.py
new file mode 100644
index 0000000..02c9183
--- /dev/null
+++ b/source/bind/python/tests/test_shock_last_valid.py
@@ -0,0 +1,67 @@
+import re
+import warnings
+
+import numpy as np
+import pytest
+
+
+_LAST_VALID_U1_CANDIDATES = (
+    1200.0,
+    1300.0,
+    1400.0,
+    1500.0,
+    1600.0,
+    1700.0,
+    1800.0,
+    1900.0,
+    2000.0,
+    2100.0,
+    2200.0,
+    2300.0,
+)
+
+
+def _find_last_valid_incident_case(cea):
+    reac = cea.Mixture(["N2 ", "CH4"], ions=True)
+    prod = cea.Mixture(["N2 ", "CH4"], products_from_reactants=True, ions=True)
+    solver = cea.ShockSolver(prod, reactants=reac, ions=True, transport=True, trace=1.0e-10)
+
+    weights = reac.moles_to_weights(np.array([0.985, 0.015]))
+    p0 = cea.units.mmhg_to_bar(0.224)
+    T0 = 291.0
+    outcomes = []
+
+    for u1 in _LAST_VALID_U1_CANDIDATES:
+        soln = cea.ShockSolution(solver, reflected=False)
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always")
+            solver.solve(soln, weights, T0, p0, u1=u1, reflected=False)
+        outcomes.append((u1, soln.last_error, soln.converged))
+        if soln.last_error == cea.LAST_VALID_SOLUTION:
+            return soln, caught
+
+    details = ", ".join(
+        f"u1={u1:.0f}: last_error={last_error}, converged={converged}"
+        for u1, last_error, converged in outcomes
+    )
+    pytest.fail(f"No LAST_VALID_SOLUTION case found across candidate velocities. {details}")
+
+
+def test_shock_last_valid_solution_warns_and_preserves_state(cea_module):
+    cea = cea_module
+    soln, caught = _find_last_valid_incident_case(cea)
+
+    pattern = re.compile(
+        r"CEA_LAST_VALID_SOLUTION.*last valid shock state retained.*converged is False"
+    )
+    assert any(
+        isinstance(w.message, RuntimeWarning) and pattern.search(str(w.message))
+        for w in caught
+    )
+
+    assert soln.last_error == cea.LAST_VALID_SOLUTION
+    assert not soln.converged
+    assert soln.P[1] > 0.0
+    assert soln.T[1] > 0.0
+    assert soln.velocity[1] > 0.0
+    assert soln.sonic_velocity[1] > 0.0
diff --git a/source/shock.f90 b/source/shock.f90
index 145a732..00e07d5 100644
--- a/source/shock.f90
+++ b/source/shock.f90
@@ -9,6 +9,10 @@ module cea_shock
     use fb_utils
     implicit none
 
+    integer, parameter :: SHOCK_SOLVE_STATUS_NOT_CONVERGED = 0
+    integer, parameter :: SHOCK_SOLVE_STATUS_CONVERGED = 1
+    integer, parameter :: SHOCK_SOLVE_STATUS_LAST_VALID = 2
+
     type :: ShockSolver
         !! Shock solver object
 
@@ -68,6 +72,9 @@ module cea_shock
         ! Convergence variables
         logical :: converged = .false.
             !! Convergence flag
+        integer :: solve_status = SHOCK_SOLVE_STATUS_NOT_CONVERGED
+            !! Internal solve outcome used by bindings to distinguish strict
+            !! non-convergence from a retained last-valid incident state.
 
     end type
     interface ShockSolution
@@ -140,6 +147,7 @@ subroutine ShockSolver_update_solution(self, soln, X1, X2, p21, t21, iter)
             t21 = EXP(log(t21) + X2)
         else  ! Converged
             soln%converged = .true.
+            soln%solve_status = SHOCK_SOLVE_STATUS_CONVERGED
         end if
 
     end subroutine
@@ -167,6 +175,7 @@ subroutine ShockSolver_fail_state(soln, idx, message)
         call log_warning("Shock calculation stopped after the last valid point.")
 
         soln%converged = .false.
+        soln%solve_status = SHOCK_SOLVE_STATUS_NOT_CONVERGED
         soln%pressure(idx) = 0.0d0
         soln%mach(idx) = 0.0d0
         soln%u(idx) = 0.0d0
@@ -405,6 +414,7 @@ subroutine ShockSolver_solve_incident(self, soln, weights, T0, P0)
         idx = 2
         G = 0.d0
         soln%converged = .false.
+        soln%solve_status = SHOCK_SOLVE_STATUS_NOT_CONVERGED
         have_last_valid = .false.
         u1 = soln%u(1)
         mach1 = soln%mach(1)
@@ -537,6 +547,11 @@ subroutine ShockSolver_solve_incident(self, soln, weights, T0, P0)
         ! singular-recovery exits and continue with warning semantics.
         if (.not. soln%converged) then
             if (ShockSolver_state_valid(soln, idx)) then
+                ! Preserve the retained incident state for inspection, but do
+                ! not mark it as converged; callers must check converged before
+                ! treating this point as a successful shock solve.
+                soln%solve_status = SHOCK_SOLVE_STATUS_LAST_VALID
+                soln%pressure(idx) = p21*P0
                 soln%rho12 = rho12
                 soln%p21 = p21
                 soln%t21 = t21
@@ -584,6 +599,8 @@ subroutine ShockSolver_solve_incident_frozen(self, soln, weights, T0, P0)
         ! Initialize
         idx = 2
         G = 0.d0
+        soln%converged = .false.
+        soln%solve_status = SHOCK_SOLVE_STATUS_NOT_CONVERGED
         u1 = soln%u(1)
         mach1 = soln%mach(1)
         soln%eq_soln(idx) = EqSolution(self%eq_solver, T_init=T0)
@@ -736,6 +753,7 @@ subroutine ShockSolver_solve_reflected(self, soln, weights, T0, P0)
         idx = 3
         G = 0.0d0  ! Reset the matrix
         soln%converged = .false.
+        soln%solve_status = SHOCK_SOLVE_STATUS_NOT_CONVERGED
         soln%eq_soln(idx) = EqSolution(self%eq_solver, T_init=soln%eq_soln(2)%T, nj_init=soln%eq_soln(2)%nj)
 
         ! Retrieve values from the incident solution
@@ -875,6 +893,7 @@ subroutine ShockSolver_solve_reflected_frozen(self, soln, weights, T0, P0)
         soln%eq_partials(idx)%dlnV_dlnP = -1.0d0
         soln%eq_partials(idx)%dlnV_dlnT = 1.0d0
         soln%converged = .false.
+        soln%solve_status = SHOCK_SOLVE_STATUS_NOT_CONVERGED
 
         ! Set the reactant weights as the species amount
         soln%eq_soln(idx)%converged = .true.
@@ -1056,6 +1075,7 @@ function ShockSolver_solve(self, reactant_weights, T0, P0, u1, mach1, reflected,
 
         ! Initialize the solution
         soln = ShockSolution_init(npts)
+        soln%solve_status = SHOCK_SOLVE_STATUS_NOT_CONVERGED
 
         ! Solve the problem
         ! --------------------------------------------------------------------
@@ -1152,6 +1172,7 @@ function ShockSolution_init(num_pts) result(self)
         integer, intent(in) :: num_pts
 
         self%num_pts = num_pts
+        self%solve_status = SHOCK_SOLVE_STATUS_NOT_CONVERGED
 
         ! Allocate data structures
         allocate(self%eq_soln(num_pts))
diff --git a/source/shock_test.pf b/source/shock_test.pf
index 96164b9..d4e221f 100644
--- a/source/shock_test.pf
+++ b/source/shock_test.pf
@@ -249,4 +249,49 @@ contains
 
     end subroutine
 
+    @test
+    subroutine test_shock_eq_last_valid_incident_status
+
+        type(Mixture) :: products
+        type(Mixture) :: reactants
+        type(ShockSolver) :: solver
+        type(ShockSolution) :: soln
+        character(:), allocatable :: product_names(:)
+        real(dp) :: p_reac, t_reac, weights(2)
+        real(dp), parameter :: candidate_u1(*) = [ &
+            1200.d0, 1300.d0, 1400.d0, 1500.d0, 1600.d0, 1700.d0, &
+            1800.d0, 1900.d0, 2000.d0, 2100.d0, 2200.d0, 2300.d0 ]
+        integer :: i
+        logical :: found_last_valid
+
+        reactants = Mixture(all_thermo, ['N2 ', 'CH4'], ions=.true.)
+        product_names = reactants%get_products(all_thermo)
+        products = Mixture(all_thermo, product_names, ions=.true.)
+
+        weights = reactants%weights_from_moles([0.985d0, 0.015d0])
+        p_reac = mmhg_to_bar(0.224d0)
+        t_reac = 291.0d0
+
+        solver = ShockSolver(products, reactants, trace=1.d-10, ions=.true., all_transport=all_transport)
+        found_last_valid = .false.
+        do i = 1, size(candidate_u1)
+            soln = solver%solve(weights, t_reac, p_reac, u1=candidate_u1(i), reflected=.false., &
+                                incident_frozen=.false., reflected_frozen=.false.)
+            if ((.not. soln%converged) .and. &
+                soln%solve_status == SHOCK_SOLVE_STATUS_LAST_VALID) then
+                found_last_valid = .true.
+                exit
+            end if
+        end do
+
+        @assertTrue(found_last_valid)
+        @assertFalse(soln%converged)
+        @assertEqual(SHOCK_SOLVE_STATUS_LAST_VALID, soln%solve_status)
+        @assertTrue(soln%p21 > 0.0d0)
+        @assertTrue(soln%eq_soln(2)%T > 0.0d0)
+        @assertTrue(soln%u(2) > 0.0d0)
+        @assertTrue(soln%v_sonic(2) > 0.0d0)
+
+    end subroutine
+
 end module