diff --git a/matlab/+occbin/solve_one_constraint.m b/matlab/+occbin/solve_one_constraint.m
index 9c60f732e4ab44fbb7316139c05662601e13c9d9..b8a16fb4a11b17714cf949e7b5e397da3b756f4a 100644
--- a/matlab/+occbin/solve_one_constraint.m
+++ b/matlab/+occbin/solve_one_constraint.m
@@ -219,7 +219,7 @@ for shock_period = 1:n_shocks_periods
             % if max_check_ahead_periods<inf, enforce unconstrained regime for period larger than max_check_ahead_periods 
             binding_constraint_new=[binding.constraint_1(1:end_periods);last_indicator];
             relaxed_constraint_new = [relax.constraint_1(1:end_periods);not(last_indicator)];
-            number_of_periods_with_violations(iter) = length(find(binding_indicator -((binding_indicator | binding_constraint_new) & ~(binding_indicator & relaxed_constraint_new))));
+            number_of_periods_with_violations(iter) = sum(binding_indicator -((binding_indicator | binding_constraint_new) & ~(binding_indicator & relaxed_constraint_new)));
             regime_exit_period(iter) = max(regime_history(shock_period).regimestart);
 
             if curb_retrench   % apply Gauss-Seidel idea of slowing down the change in the guess
@@ -259,7 +259,7 @@ for shock_period = 1:n_shocks_periods
                     if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
                         %                     vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
                         %                     is_periodic(kiter) = isequal(vvv, binding_indicator);
-                        is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && length(find(binding_indicator_history{iter}-binding_indicator))<=opts_simul_.periodic_solution_threshold;
+                        is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && sum(binding_indicator_history{iter}-binding_indicator)<=opts_simul_.periodic_solution_threshold;
                     else
                         is_periodic(kiter)=false;
                     end
diff --git a/matlab/+occbin/solve_two_constraints.m b/matlab/+occbin/solve_two_constraints.m
index ba4bdfeec7fbcf6bf843fe2a80fa6687b02bf5fa..f4c24f14a5f7f712400bf12efd1d0dca894b3738 100644
--- a/matlab/+occbin/solve_two_constraints.m
+++ b/matlab/+occbin/solve_two_constraints.m
@@ -243,8 +243,8 @@ for shock_period = 1:n_shocks_periods
             
             binding_constraint_new=[binding.constraint_1(1:end_periods);last_indicator; binding.constraint_2(1:end_periods);last_indicator];
             relaxed_constraint_new = [relax.constraint_1(1:end_periods);not(last_indicator); relax.constraint_2(1:end_periods);not(last_indicator)];
-            tmp_nper(1) = length(find(binding_indicator(:,1) - (binding_indicator(:,1) | [binding.constraint_1(1:end_periods);last_indicator]) & ~(binding_indicator(:,1) & [relax.constraint_1(1:end_periods);not(last_indicator)])));
-            tmp_nper(2) = length(find(binding_indicator(:,2) - (binding_indicator(:,2) | [binding.constraint_2(1:end_periods);last_indicator]) & ~(binding_indicator(:,2) & [relax.constraint_2(1:end_periods);not(last_indicator)])));
+            tmp_nper(1) = sum(binding_indicator(:,1) - (binding_indicator(:,1) | [binding.constraint_1(1:end_periods);last_indicator]) & ~(binding_indicator(:,1) & [relax.constraint_1(1:end_periods);not(last_indicator)]));
+            tmp_nper(2) = sum(binding_indicator(:,2) - (binding_indicator(:,2) | [binding.constraint_2(1:end_periods);last_indicator]) & ~(binding_indicator(:,2) & [relax.constraint_2(1:end_periods);not(last_indicator)]));
             number_of_periods_with_violations(iter) = max(tmp_nper);
             regime_exit_period(iter,1) = max(regime_history(shock_period).regimestart1);
             regime_exit_period(iter,2) = max(regime_history(shock_period).regimestart2);
@@ -290,7 +290,7 @@ for shock_period = 1:n_shocks_periods
                     if size(binding_indicator,1)== size(binding_indicator_history{kiter},1)
                         %                     vvv = [binding_indicator_history{kiter}; false(size(binding_indicator,1)- size(binding_indicator_history{kiter},1), 1)];
                         %                     is_periodic(kiter) = isequal(vvv, binding_indicator);
-                        is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && length(find(binding_indicator_history{iter}(:,1)-binding_indicator(:,1)))<=opts_simul_.periodic_solution_threshold  && length(find(binding_indicator_history{iter}(:,2)-binding_indicator(:,2)))<=opts_simul_.periodic_solution_threshold;
+                        is_periodic(kiter) = isequal(binding_indicator_history{kiter}, binding_indicator) && sum(binding_indicator_history{iter}(:,1)-binding_indicator(:,1))<=opts_simul_.periodic_solution_threshold  && sum(binding_indicator_history{iter}(:,2)-binding_indicator(:,2))<=opts_simul_.periodic_solution_threshold;
                     else
                         is_periodic(kiter)=false;
                     end