diff --git a/R/engine_functions.R b/R/engine_functions.R
index ac536313..dc3bf9cb 100644
--- a/R/engine_functions.R
+++ b/R/engine_functions.R
@@ -874,6 +874,28 @@
 #' engine_eval(~ check_finite(1/0)) ## returns nothing and throws an error
 #' ```
 #' 
+#' ## Control Flow
+#' 
+#' Functions for controlling the flow of simulations.
+#' 
+#' ### Functions
+#' 
+#' * `stop_if_lt(x, y)` : Stop the simulation loop and return
+#' simulations if any element of `x` is less than the
+#' corresponding element in `y`. The elements of `y` are
+#' recycled to match those in `x` if necessary and if possible.
+#' 
+#' ### Arguments
+#' 
+#' * `x` : Matrix to compare with `y`.
+#' * `y` : Matrix to compare with `x`, with a default value of
+#' `0`.
+#'
+#' ### Return 
+#'
+#' An empty matrix. This function is returned for its effect of
+#' breaking out of the simulation loop.
+#'
 #' ## Assign (deprecated)
 #'
 #' Assign values to a subset of the elements in a matrix.
@@ -1056,6 +1078,7 @@
 #' @aliases pnorm
 #' @aliases invlogit
 #' @aliases logit
+#' @aliases stop_if_lt
 #' @aliases cumsum
 #' @aliases assign
 #' @aliases unpack
diff --git a/R/enum.R b/R/enum.R
index f655bb28..771d16d9 100644
--- a/R/enum.R
+++ b/R/enum.R
@@ -59,6 +59,7 @@ valid_func_sigs = c(
       , "fwrap: pnorm(q, mean, sd)"
       , "fwrap: invlogit(x)"
       , "fwrap: logit(x)"
+      , "fwrap: stop_if_lt(x, y)"
       , "fwrap: cumsum(x)"
       , "fwrap: assign(x, i, j, v)"
       , "fwrap: unpack(x, ...)"
diff --git a/R/mp_tmb_model_spec.R b/R/mp_tmb_model_spec.R
index 280becdf..5d667443 100644
--- a/R/mp_tmb_model_spec.R
+++ b/R/mp_tmb_model_spec.R
@@ -42,6 +42,7 @@ TMBModelSpec = function(
         self$update_method$before()
       , self$update_method$during()
       , self$update_method$after()
+      , .simulate_exprs = self$sim_exprs
     )
   }
   
diff --git a/inst/starter_models/nfds/README.md b/inst/starter_models/nfds/README.md
index e7e221ea..8030dba5 100644
--- a/inst/starter_models/nfds/README.md
+++ b/inst/starter_models/nfds/README.md
@@ -355,14 +355,14 @@ start_time <- Sys.time()
 old_sim = mp_simulator(old_spec,time_steps = 1L,outputs = "Y")
 end_time <- Sys.time()
 end_time - start_time
-#> Time difference of 0.2173371 secs
+#> Time difference of 0.195493 secs
 
 # new spec
 start_time <- Sys.time()
 new_sim = mp_simulator(new_spec,time_steps = 1L,outputs = "Y")
 end_time <- Sys.time()
 end_time - start_time
-#> Time difference of 0.222172 secs
+#> Time difference of 0.187727 secs
 ```
 
 Creating a simulator with one time step results in the new model
diff --git a/man/engine_functions.Rd b/man/engine_functions.Rd
index c54633a8..207aa8f2 100644
--- a/man/engine_functions.Rd
+++ b/man/engine_functions.Rd
@@ -61,6 +61,7 @@
 \alias{pnorm}
 \alias{invlogit}
 \alias{logit}
+\alias{stop_if_lt}
 \alias{cumsum}
 \alias{assign}
 \alias{unpack}
@@ -1070,6 +1071,34 @@ engine_eval(~ check_finite(1/0)) ## returns nothing and throws an error
 
 }
 
+\subsection{Control Flow}{
+
+Functions for controlling the flow of simulations.
+\subsection{Functions}{
+\itemize{
+\item \code{stop_if_lt(x, y)} : Stop the simulation loop and return
+simulations if any element of \code{x} is less than the
+corresponding element in \code{y}. The elements of \code{y} are
+recycled to match those in \code{x} if necessary and if possible.
+}
+}
+
+\subsection{Arguments}{
+\itemize{
+\item \code{x} : Matrix to compare with \code{y}.
+\item \code{y} : Matrix to compare with \code{x}, with a default value of
+\code{0}.
+}
+}
+
+\subsection{Return}{
+
+An empty matrix. This function is returned for its effect of
+breaking out of the simulation loop.
+}
+
+}
+
 \subsection{Assign (deprecated)}{
 
 Assign values to a subset of the elements in a matrix.
diff --git a/misc/dev/dev.cpp b/misc/dev/dev.cpp
index 179d1d4a..b2a5b04e 100644
--- a/misc/dev/dev.cpp
+++ b/misc/dev/dev.cpp
@@ -135,9 +135,10 @@ enum macpan2_func
     , MP2_PNORM = 57 // fwrap: pnorm(q, mean, sd)
     , MP2_INVLOGIT = 58 // fwrap: invlogit(x)
     , MP2_LOGIT = 59 // fwrap: logit(x)
-    , MP2_CUMSUM = 60 // fwrap: cumsum(x)
-    , MP2_ASSIGN = 61 // fwrap: assign(x, i, j, v)
-    , MP2_UNPACK = 62 // fwrap: unpack(x, ...)
+    , MP2_STOP_IF_LT = 60 // fwrap: stop_if_lt(x, y)
+    , MP2_CUMSUM = 61 // fwrap: cumsum(x)
+    , MP2_ASSIGN = 62 // fwrap: assign(x, i, j, v)
+    , MP2_UNPACK = 63 // fwrap: unpack(x, ...)
 };
 
 enum macpan2_meth
@@ -825,7 +826,6 @@ class ArgList
                     // Rcpp::Rcout << "step f" << std::endl;
                     // Rcpp::Rcout << "bad column vector" << std::endl;
                     error_code = 501;
-                    // break; // Exit the loop on error
                 }
             }
             else if (mat.cols() == cols) {
@@ -839,7 +839,6 @@ class ArgList
                 {
                     // Rcpp::Rcout << "bad row vector" << std::endl;
                     error_code = 501;
-                    // break; // Exit the loop on error
                 }
             }
             else
@@ -847,7 +846,6 @@ class ArgList
                 // Rcpp::Rcout << "really bad" << std::endl;
                 // Rcpp::Rcout << "step g" << std::endl;
                 error_code = 501;
-                // break; // Exit the loop on error
             }
 
             if (error_code != 0)
@@ -901,11 +899,6 @@ class ArgList
     int size_;
     int error_code_ = 0; // Initialize the error code to 0 (no error) by default
 };
-// private:
-//     std::vector<ItemType> items_;
-//     int size_;
-//     int error_code_ = 0; // Initialize the error code to 0 (no error) by default
-// };
 
 template <class Type>
 class ExprEvaluator
@@ -920,6 +913,7 @@ class ExprEvaluator
     ListOfIntVecs meth_int_vecs;
     ListOfIntVecs valid_int_vecs;
     vector<Type> valid_literals;
+    bool early_stop;
 
 public:
     // constructor
@@ -951,6 +945,10 @@ class ExprEvaluator
         meth_int_vecs = meth_int_vecs_;
         valid_int_vecs = valid_int_vecs_;
         valid_literals = valid_literals_;
+        // Rcpp::Rcout << "setting early_stop to false" << std::endl;
+        early_stop = false;
+        // Rcpp::Rcout << "early_stop now equals: " << early_stop << std::endl;
+        // Rcpp::Rcout << "but GetEarlyStop() now returns: " << GetEarlyStop() << std::endl;
 
         strcpy(error_message, "OK");
     };
@@ -964,6 +962,7 @@ class ExprEvaluator
     std::vector<int> GetArgRows() { return arg_rows; };
     std::vector<int> GetArgCols() { return arg_cols; };
     std::vector<int> GetArgTypeInts() { return arg_type_ints; };
+    bool GetEarlyStop() { return early_stop; };
 
     // setters
     void SetError(unsigned char code, const char *message, int e_row, int f_int, std::vector<int> a_rows, std::vector<int> a_cols, std::vector<int> a_type_ints, int t_int) {
@@ -976,13 +975,17 @@ class ExprEvaluator
         arg_type_ints = a_type_ints;
         strcpy(error_message, message);
     };
+    void SetEarlyStop() {
+        early_stop = true;
+    };
 
     // evaluator
     matrix<Type> EvalExpr(
         const vector<ListOfMatrices<Type>> &hist, // current simulation history
         int t,                                    // current time step
         ListOfMatrices<Type> &valid_vars,         // current list of values of each matrix
-        int row = 0                               // current expression parse table row being evaluated
+        int row = 0,                              // current expression parse table row being evaluated
+        bool simulate_block = false               // evaluated within simulate block?
     ) {
 
         // total number of time steps in the simulation loop
@@ -1124,7 +1127,7 @@ class ExprEvaluator
                 {
                     // otherwise, recursively descend into the parse tree
                     // to pick out the arguments
-                    args.set(i, EvalExpr(hist, t, valid_vars, table_i[row] + i));
+                    args.set(i, EvalExpr(hist, t, valid_vars, table_i[row] + i, simulate_block));
                 }
 
                 // Check here if index2mats actually points at
@@ -3144,6 +3147,50 @@ class ExprEvaluator
                 }
                 return m;
             
+            
+            // #' ## Control Flow
+            // #' 
+            // #' Functions for controlling the flow of simulations.
+            // #' 
+            // #' ### Functions
+            // #' 
+            // #' * `stop_if_lt(x, y)` : Stop the simulation loop and return
+            // #' simulations if any element of `x` is less than the
+            // #' corresponding element in `y`. The elements of `y` are
+            // #' recycled to match those in `x` if necessary and if possible.
+            // #' 
+            // #' ### Arguments
+            // #' 
+            // #' * `x` : Matrix to compare with `y`.
+            // #' * `y` : Matrix to compare with `x`, with a default value of
+            // #' `0`.
+            // #'
+            // #' ### Return 
+            // #'
+            // #' An empty matrix. This function is returned for its effect of
+            // #' breaking out of the simulation loop.
+            // #'
+            case MP2_STOP_IF_LT:
+                m = args.get_as_mat(0);
+                rows = m.rows();
+                cols = m.cols();
+                if (n == 1) {
+                    m1 = matrix<Type>::Zero(rows, cols);
+                } else {
+                    v1.push_back(1);
+                    args = args.recycle_to_shape(v1, rows, cols);
+                    m1 = args.get_as_mat(1);
+                }
+                for (int i = 0; i < rows; i++) {
+                    for (int j = 0; j < cols; j++) {
+                        if (m.coeff(i, j) < m1.coeff(i, j)) {
+                            SetEarlyStop();
+                            return m2; // empty matrix
+                        }
+                    }
+                }
+                return m2; // empty matrix
+            
             // #' ## Assign (deprecated)
             // #'
             // #' Assign values to a subset of the elements in a matrix.
@@ -3250,7 +3297,7 @@ class ExprEvaluator
                     valid_vars.m_matrices[matIndex].coeffRef(rowIndex, colIndex) = args[3].coeff(k, 0);
                 }
                 return m2; // empty matrix
-
+            
             // #' ## Unpack (deprecated)
             // #'
             // #' Unpack elements of a matrix into smaller matrices.
@@ -3529,6 +3576,7 @@ class MatAssigner {
         int func_int = exprEvaluator.GetFuncInt();               \
         int time_int = exprEvaluator.GetTimeInt();               \
         const char *err_msg = exprEvaluator.GetErrorMessage();   \
+        bool early_stop = exprEvaluator.GetEarlyStop();          \
         REPORT(error);                                           \
         REPORT(expr_row);                                        \
         REPORT(func_int);                                        \
@@ -3543,6 +3591,7 @@ class MatAssigner {
         logfile << "Expression row = " << expr_row << std::endl; \
         logfile << "Function code = " << func_int << std::endl;  \
         logfile << "Time step = " << time_int << std::endl;      \
+        logfile << "Early stop = " << early_stop << std::endl;   \
         logfile.close();                                         \
     }
 
@@ -3719,6 +3768,14 @@ Type objective_function<Type>::operator()()
     for (int i = 0; i < n; i++)
         mats.m_matrices[r_mat_id[i]].coeffRef(r_row_id[i], r_col_id[i]) = random[r_par_id[i]];
 
+    // in case you need to stop early you can use this to fill the remaining 
+    // iterations of the simulation history
+    ListOfMatrices<Type> empty_mats_list(mats);
+    matrix<Type> empty_matrix;
+    for (int i = 0; i < empty_mats_list.m_matrices.size(); i++) {
+        empty_mats_list.m_matrices[i] = empty_matrix;
+    }
+    
     // Simulation history
     /// each element of this history 'vector' is a list of the matrices
     /// in the model at a particular point in history
@@ -3773,7 +3830,7 @@ Type objective_function<Type>::operator()()
         if (expr_sim_block[i] == 1) {
             SIMULATE {
                 result = exprEvaluator.EvalExpr(
-                    simulation_history, 0, mats, p_table_row);
+                    simulation_history, 0, mats, p_table_row, true);
             }
         }
         else
@@ -3821,23 +3878,34 @@ Type objective_function<Type>::operator()()
             Rcpp::Rcout << "Eval expression --- " << i << std::endl;
             Rcpp::Rcout << "expr_num_p_table_rows[i] " << expr_num_p_table_rows[expr_index + i] << std::endl;
 #endif
-            matrix<Type> result;
-            if (expr_sim_block[i] == 1) {
-                SIMULATE {
+            if (!exprEvaluator.GetEarlyStop()) {
+                matrix<Type> result;
+                if (expr_sim_block[i] == 1) {
+                    SIMULATE {
+                        result = exprEvaluator.EvalExpr(
+                              simulation_history
+                            , k + 1
+                            , mats
+                            , p_table_row2
+                            , true
+                        );
+                    }
+                }
+                else
                     result = exprEvaluator.EvalExpr(
-                        simulation_history, k + 1, mats, p_table_row2);
+                          simulation_history
+                        , k + 1
+                        , mats
+                        , p_table_row2
+                        , false
+                    );
+    
+                if (exprEvaluator.GetErrorCode()) {
+                    REPORT_ERROR
+                    return 0.0;
                 }
+                matAssigner.matAssign(result, mats, a_table_row2);
             }
-            else
-                result = exprEvaluator.EvalExpr(
-                    simulation_history, k + 1, mats, p_table_row2);
-
-            if (exprEvaluator.GetErrorCode()) {
-                REPORT_ERROR
-                return 0.0;
-            }
-            // mats.m_matrices[expr_output_id[expr_index+i]] = result;
-            matAssigner.matAssign(result, mats, a_table_row2);
 
             p_table_row2 += expr_num_p_table_rows[expr_index + i];
             a_table_row2 += assign_num_a_table_rows[expr_index + i];
@@ -3848,13 +3916,24 @@ Type objective_function<Type>::operator()()
                 Rcpp::Rcout << "mats = " << mats.m_matrices[ii] << std::endl;
 #endif
         }
-        // simulation_history[k+1] = mats;
-        UpdateSimulationHistory(
-            simulation_history,
-            k + 1,
-            mats,
-            mats_save_hist,
-            hist_shape_template);
+        if (exprEvaluator.GetEarlyStop()) {
+            UpdateSimulationHistory(
+                simulation_history,
+                k + 1,
+                empty_mats_list,
+                mats_save_hist,
+                hist_shape_template
+            );
+        } else {
+            // simulation_history[k+1] = mats;
+            UpdateSimulationHistory(
+                simulation_history,
+                k + 1,
+                mats,
+                mats_save_hist,
+                hist_shape_template
+            );
+        }
     }
     p_table_row = p_table_row2;
     a_table_row = a_table_row2;
@@ -3871,7 +3950,7 @@ Type objective_function<Type>::operator()()
         if (expr_sim_block[i] == 1) {
             SIMULATE {
                 result = exprEvaluator.EvalExpr(
-                    simulation_history, time_steps + 1, mats, p_table_row);
+                    simulation_history, time_steps + 1, mats, p_table_row, true);
             }
         }
         else {
@@ -3898,6 +3977,7 @@ Type objective_function<Type>::operator()()
         mats,
         mats_save_hist,
         hist_shape_template);
+    
 
 #ifdef MP_VERBOSE
     Rcpp::Rcout << "Simulation history ..." << std::endl;
@@ -3969,6 +4049,7 @@ Type objective_function<Type>::operator()()
         matrix<Type> value_column = values.block(0, 4, values.rows(), 1);
         ADREPORT(value_column)
     }
+    
 
     // 7 Calc the return of the objective function
     matrix<Type> ret;
@@ -3985,5 +4066,6 @@ Type objective_function<Type>::operator()()
 #ifdef MP_VERBOSE
     Rcpp::Rcout << "======== end of objective function ========" << std::endl;
 #endif
+    
     return ret.coeff(0, 0);
 }
diff --git a/src/macpan2.cpp b/src/macpan2.cpp
index 54ae64d5..3a55c426 100644
--- a/src/macpan2.cpp
+++ b/src/macpan2.cpp
@@ -136,9 +136,10 @@ enum macpan2_func
     , MP2_PNORM = 57 // fwrap: pnorm(q, mean, sd)
     , MP2_INVLOGIT = 58 // fwrap: invlogit(x)
     , MP2_LOGIT = 59 // fwrap: logit(x)
-    , MP2_CUMSUM = 60 // fwrap: cumsum(x)
-    , MP2_ASSIGN = 61 // fwrap: assign(x, i, j, v)
-    , MP2_UNPACK = 62 // fwrap: unpack(x, ...)
+    , MP2_STOP_IF_LT = 60 // fwrap: stop_if_lt(x, y)
+    , MP2_CUMSUM = 61 // fwrap: cumsum(x)
+    , MP2_ASSIGN = 62 // fwrap: assign(x, i, j, v)
+    , MP2_UNPACK = 63 // fwrap: unpack(x, ...)
 };
 
 enum macpan2_meth
@@ -826,7 +827,6 @@ class ArgList
                     // Rcpp::Rcout << "step f" << std::endl;
                     // Rcpp::Rcout << "bad column vector" << std::endl;
                     error_code = 501;
-                    // break; // Exit the loop on error
                 }
             }
             else if (mat.cols() == cols) {
@@ -840,7 +840,6 @@ class ArgList
                 {
                     // Rcpp::Rcout << "bad row vector" << std::endl;
                     error_code = 501;
-                    // break; // Exit the loop on error
                 }
             }
             else
@@ -848,7 +847,6 @@ class ArgList
                 // Rcpp::Rcout << "really bad" << std::endl;
                 // Rcpp::Rcout << "step g" << std::endl;
                 error_code = 501;
-                // break; // Exit the loop on error
             }
 
             if (error_code != 0)
@@ -902,11 +900,6 @@ class ArgList
     int size_;
     int error_code_ = 0; // Initialize the error code to 0 (no error) by default
 };
-// private:
-//     std::vector<ItemType> items_;
-//     int size_;
-//     int error_code_ = 0; // Initialize the error code to 0 (no error) by default
-// };
 
 template <class Type>
 class ExprEvaluator
@@ -921,6 +914,7 @@ class ExprEvaluator
     ListOfIntVecs meth_int_vecs;
     ListOfIntVecs valid_int_vecs;
     vector<Type> valid_literals;
+    bool early_stop;
 
 public:
     // constructor
@@ -952,6 +946,10 @@ class ExprEvaluator
         meth_int_vecs = meth_int_vecs_;
         valid_int_vecs = valid_int_vecs_;
         valid_literals = valid_literals_;
+        // Rcpp::Rcout << "setting early_stop to false" << std::endl;
+        early_stop = false;
+        // Rcpp::Rcout << "early_stop now equals: " << early_stop << std::endl;
+        // Rcpp::Rcout << "but GetEarlyStop() now returns: " << GetEarlyStop() << std::endl;
 
         strcpy(error_message, "OK");
     };
@@ -965,6 +963,7 @@ class ExprEvaluator
     std::vector<int> GetArgRows() { return arg_rows; };
     std::vector<int> GetArgCols() { return arg_cols; };
     std::vector<int> GetArgTypeInts() { return arg_type_ints; };
+    bool GetEarlyStop() { return early_stop; };
 
     // setters
     void SetError(unsigned char code, const char *message, int e_row, int f_int, std::vector<int> a_rows, std::vector<int> a_cols, std::vector<int> a_type_ints, int t_int) {
@@ -977,13 +976,17 @@ class ExprEvaluator
         arg_type_ints = a_type_ints;
         strcpy(error_message, message);
     };
+    void SetEarlyStop() {
+        early_stop = true;
+    };
 
     // evaluator
     matrix<Type> EvalExpr(
         const vector<ListOfMatrices<Type>> &hist, // current simulation history
         int t,                                    // current time step
         ListOfMatrices<Type> &valid_vars,         // current list of values of each matrix
-        int row = 0                               // current expression parse table row being evaluated
+        int row = 0,                              // current expression parse table row being evaluated
+        bool simulate_block = false               // evaluated within simulate block?
     ) {
 
         // total number of time steps in the simulation loop
@@ -1125,7 +1128,7 @@ class ExprEvaluator
                 {
                     // otherwise, recursively descend into the parse tree
                     // to pick out the arguments
-                    args.set(i, EvalExpr(hist, t, valid_vars, table_i[row] + i));
+                    args.set(i, EvalExpr(hist, t, valid_vars, table_i[row] + i, simulate_block));
                 }
 
                 // Check here if index2mats actually points at
@@ -3145,6 +3148,50 @@ class ExprEvaluator
                 }
                 return m;
             
+            
+            // #' ## Control Flow
+            // #' 
+            // #' Functions for controlling the flow of simulations.
+            // #' 
+            // #' ### Functions
+            // #' 
+            // #' * `stop_if_lt(x, y)` : Stop the simulation loop and return
+            // #' simulations if any element of `x` is less than the
+            // #' corresponding element in `y`. The elements of `y` are
+            // #' recycled to match those in `x` if necessary and if possible.
+            // #' 
+            // #' ### Arguments
+            // #' 
+            // #' * `x` : Matrix to compare with `y`.
+            // #' * `y` : Matrix to compare with `x`, with a default value of
+            // #' `0`.
+            // #'
+            // #' ### Return 
+            // #'
+            // #' An empty matrix. This function is returned for its effect of
+            // #' breaking out of the simulation loop.
+            // #'
+            case MP2_STOP_IF_LT:
+                m = args.get_as_mat(0);
+                rows = m.rows();
+                cols = m.cols();
+                if (n == 1) {
+                    m1 = matrix<Type>::Zero(rows, cols);
+                } else {
+                    v1.push_back(1);
+                    args = args.recycle_to_shape(v1, rows, cols);
+                    m1 = args.get_as_mat(1);
+                }
+                for (int i = 0; i < rows; i++) {
+                    for (int j = 0; j < cols; j++) {
+                        if (m.coeff(i, j) < m1.coeff(i, j)) {
+                            SetEarlyStop();
+                            return m2; // empty matrix
+                        }
+                    }
+                }
+                return m2; // empty matrix
+            
             // #' ## Assign (deprecated)
             // #'
             // #' Assign values to a subset of the elements in a matrix.
@@ -3251,7 +3298,7 @@ class ExprEvaluator
                     valid_vars.m_matrices[matIndex].coeffRef(rowIndex, colIndex) = args[3].coeff(k, 0);
                 }
                 return m2; // empty matrix
-
+            
             // #' ## Unpack (deprecated)
             // #'
             // #' Unpack elements of a matrix into smaller matrices.
@@ -3530,6 +3577,7 @@ class MatAssigner {
         int func_int = exprEvaluator.GetFuncInt();               \
         int time_int = exprEvaluator.GetTimeInt();               \
         const char *err_msg = exprEvaluator.GetErrorMessage();   \
+        bool early_stop = exprEvaluator.GetEarlyStop();          \
         REPORT(error);                                           \
         REPORT(expr_row);                                        \
         REPORT(func_int);                                        \
@@ -3544,6 +3592,7 @@ class MatAssigner {
         logfile << "Expression row = " << expr_row << std::endl; \
         logfile << "Function code = " << func_int << std::endl;  \
         logfile << "Time step = " << time_int << std::endl;      \
+        logfile << "Early stop = " << early_stop << std::endl;   \
         logfile.close();                                         \
     }
 
@@ -3720,6 +3769,14 @@ Type objective_function<Type>::operator()()
     for (int i = 0; i < n; i++)
         mats.m_matrices[r_mat_id[i]].coeffRef(r_row_id[i], r_col_id[i]) = random[r_par_id[i]];
 
+    // in case you need to stop early you can use this to fill the remaining 
+    // iterations of the simulation history
+    ListOfMatrices<Type> empty_mats_list(mats);
+    matrix<Type> empty_matrix;
+    for (int i = 0; i < empty_mats_list.m_matrices.size(); i++) {
+        empty_mats_list.m_matrices[i] = empty_matrix;
+    }
+    
     // Simulation history
     /// each element of this history 'vector' is a list of the matrices
     /// in the model at a particular point in history
@@ -3774,7 +3831,7 @@ Type objective_function<Type>::operator()()
         if (expr_sim_block[i] == 1) {
             SIMULATE {
                 result = exprEvaluator.EvalExpr(
-                    simulation_history, 0, mats, p_table_row);
+                    simulation_history, 0, mats, p_table_row, true);
             }
         }
         else
@@ -3822,23 +3879,34 @@ Type objective_function<Type>::operator()()
             Rcpp::Rcout << "Eval expression --- " << i << std::endl;
             Rcpp::Rcout << "expr_num_p_table_rows[i] " << expr_num_p_table_rows[expr_index + i] << std::endl;
 #endif
-            matrix<Type> result;
-            if (expr_sim_block[i] == 1) {
-                SIMULATE {
+            if (!exprEvaluator.GetEarlyStop()) {
+                matrix<Type> result;
+                if (expr_sim_block[i] == 1) {
+                    SIMULATE {
+                        result = exprEvaluator.EvalExpr(
+                              simulation_history
+                            , k + 1
+                            , mats
+                            , p_table_row2
+                            , true
+                        );
+                    }
+                }
+                else
                     result = exprEvaluator.EvalExpr(
-                        simulation_history, k + 1, mats, p_table_row2);
+                          simulation_history
+                        , k + 1
+                        , mats
+                        , p_table_row2
+                        , false
+                    );
+    
+                if (exprEvaluator.GetErrorCode()) {
+                    REPORT_ERROR
+                    return 0.0;
                 }
+                matAssigner.matAssign(result, mats, a_table_row2);
             }
-            else
-                result = exprEvaluator.EvalExpr(
-                    simulation_history, k + 1, mats, p_table_row2);
-
-            if (exprEvaluator.GetErrorCode()) {
-                REPORT_ERROR
-                return 0.0;
-            }
-            // mats.m_matrices[expr_output_id[expr_index+i]] = result;
-            matAssigner.matAssign(result, mats, a_table_row2);
 
             p_table_row2 += expr_num_p_table_rows[expr_index + i];
             a_table_row2 += assign_num_a_table_rows[expr_index + i];
@@ -3849,13 +3917,24 @@ Type objective_function<Type>::operator()()
                 Rcpp::Rcout << "mats = " << mats.m_matrices[ii] << std::endl;
 #endif
         }
-        // simulation_history[k+1] = mats;
-        UpdateSimulationHistory(
-            simulation_history,
-            k + 1,
-            mats,
-            mats_save_hist,
-            hist_shape_template);
+        if (exprEvaluator.GetEarlyStop()) {
+            UpdateSimulationHistory(
+                simulation_history,
+                k + 1,
+                empty_mats_list,
+                mats_save_hist,
+                hist_shape_template
+            );
+        } else {
+            // simulation_history[k+1] = mats;
+            UpdateSimulationHistory(
+                simulation_history,
+                k + 1,
+                mats,
+                mats_save_hist,
+                hist_shape_template
+            );
+        }
     }
     p_table_row = p_table_row2;
     a_table_row = a_table_row2;
@@ -3872,7 +3951,7 @@ Type objective_function<Type>::operator()()
         if (expr_sim_block[i] == 1) {
             SIMULATE {
                 result = exprEvaluator.EvalExpr(
-                    simulation_history, time_steps + 1, mats, p_table_row);
+                    simulation_history, time_steps + 1, mats, p_table_row, true);
             }
         }
         else {
@@ -3899,6 +3978,7 @@ Type objective_function<Type>::operator()()
         mats,
         mats_save_hist,
         hist_shape_template);
+    
 
 #ifdef MP_VERBOSE
     Rcpp::Rcout << "Simulation history ..." << std::endl;
@@ -3970,6 +4050,7 @@ Type objective_function<Type>::operator()()
         matrix<Type> value_column = values.block(0, 4, values.rows(), 1);
         ADREPORT(value_column)
     }
+    
 
     // 7 Calc the return of the objective function
     matrix<Type> ret;
@@ -3986,5 +4067,6 @@ Type objective_function<Type>::operator()()
 #ifdef MP_VERBOSE
     Rcpp::Rcout << "======== end of objective function ========" << std::endl;
 #endif
+    
     return ret.coeff(0, 0);
 }
diff --git a/tests/testthat/test-stop-if.R b/tests/testthat/test-stop-if.R
new file mode 100644
index 00000000..8e80c858
--- /dev/null
+++ b/tests/testthat/test-stop-if.R
@@ -0,0 +1,16 @@
+test_that("simulations can be stopped early", {
+  exprs = list(x ~ x - 2, sims = dummy ~ stop_if_lt(x), y ~ y - 2)
+  def = list(x = 10, y = 10)
+  ts = 10
+  spec = mp_tmb_model_spec(during = exprs, default = def)
+  sim = mp_simulator(spec, ts, c("x", "y"))
+  
+  expect_equal(
+      mp_final(sim)$value
+    , c(-2, 0)
+  )
+  expect_equal(
+      mp_trajectory(sim)$value
+    , rep(seq(8, 0, by = -2), each = 2)
+  )
+})