diff --git a/PF_insert.m b/PF_insert.m
index fa19bd8..3ef6262 100644
--- a/PF_insert.m
+++ b/PF_insert.m
@@ -1,6 +1,8 @@
 %PF_insert. Guess and verify part of the algorithm.  Written by Michael Hatcher (m.c.hatcher@soton.ac.uk). 
 %Any errors are my own. Last updated: 12/04/2023.
 
+%Housekeeping
+mstar = zeros(N_guess,1);
 X_stack = NaN(nvar,1);
 X_sol = NaN(nvar,T_sim-1,N_guess); X_sol_exc = NaN(nvar-nx,T_sim-1,N_guess);
 ind_sol = NaN(N_guess,T_sim-1); 
@@ -78,23 +80,27 @@
         X_star = F*[X; X_e; X_lag] + G*e_vec(:,t) + H;
         X_star_store(t) = X_star;
         
-        if X_star >= X1_min && ind(t) == 1  || X_star <= X1_min && ind(t) == 0
+        if X_star > X1_min && ind(t) == 1  || X_star <= X1_min && ind(t) == 0    %Can use X_star => X1_min && ind(t) == 1 for knife-edge case
             Verify(t) = 1;
         else
             break
         end
 
-        %if X_stack(1,t) == max(X1_min,X_star) does not work well due to numerical precision
-        %For accuracy, can use vpa(.)
+        %if X_stack(1,t) == max(X1_min,X_star) does not seem to work well
             
 end 
 
-    if sum(Verify) == T_sim-1 && min(X_stack(1,T_guess:T_sim-1)) > X1_min  
+    if sum(Verify) == T_sim-1 && min(X_stack(1,T_guess+1:T_sim-1)) > X1_min  
          X_sol(:,:,m) = X_stack(:,:);
          X_sol_exc(:,:,m) = X_stack(1:nvar-nx,:);
          mstar(m) = m;
          ind_sol(m,:) = ind(1:T_sim-1)';
          X_star_sol(m,:) = X_star_store(1:T_sim-1);
+
+         if ~isempty(not_P) && not_P==0  %          
+            break    %Stop search if M is a P-matrix and a solution is found
+        end
+
     end
          
 end
diff --git a/P_matrix.m b/P_matrix.m
index 1c8b007..bfa8a45 100644
--- a/P_matrix.m
+++ b/P_matrix.m
@@ -1,18 +1,19 @@
-%Script to determine if M is a P-matrix. %Uses ptest.m written by Michael
-%Tsatsomeros (https://www.math.wsu.edu/faculty/tsat/matlab.html). The code can also be used with ptest3.m, by the same author.
-%Written by Michael Hatcher (m.c.hatcher@soton.ac.uk). Any errors are my own. Last updated: 12/04/2023.
+%Script to determine if M is a P-matrix. %Uses ptest.m and ptest3.m written by Michael Tsatsomeros 
+%(see the MATLAB depot at https://www.math.wsu.edu/faculty/tsat/matlab.html). 
+%This script is written by Michael Hatcher (m.c.hatcher@soton.ac.uk). Any errors are my own. Last updated: 12/04/2023.
 
 run Pos_Def.m
 
-if not_P == 0 && is_pd == 0       
+if not_P ~=1 && is_pd == 0   
+
+    %ptest_val = ptest(M);
+    ptest_val = ptest3(M);  %faster version
     
-    if ptest(M)==0  
+    if ptest_val==0  
         not_P = 1;
+    else
+        not_P = 0;
     end
-    
-    %if ptest3(M)==0  
-    %    not_P = 1;
-    %end
 
 end
 
diff --git a/Pos_Def.m b/Pos_Def.m
index 5ff7025..2d18309 100644
--- a/Pos_Def.m
+++ b/Pos_Def.m
@@ -1,35 +1,18 @@
 %Check if M matrix is general positive definite --> P-matrix
-%Check if simple necessary conditions for a P-matrix are violated
 %Written by Michael Hatcher (m.c.hatcher@soton.ac.uk). 
-%Any errors are my own. Last updated: 12/04/2023.
+%Any errors are my own. Last updated: 03/10/2024.
 
-%Based on Holden (2022, Lemma 1 Part 2)
+%Based on Holden (2023, Lemma 1 Part 2)
 eig_vals = eig(M + M'); 
 length_M = length(M);
-is_pd = nnz(eig_vals>0) == length_M; % flag to check if it is PD 
+is_pd = nnz(eig_vals>0) == length_M; % flag to check if it is PD, nnz - no. of non-zero elements
 
-if is_pd == 0
+if is_pd == 1
 
-    if det(M) <= 0 || min(diag(M)) <=0 
-        not_P = 1;
+    not_P = 0;
 
-    else  
+elseif min(diag(M)) <=0  || det(M) <= 0 
 
-n_skip = max(round(length_M/10),1);
-max_iter = max(length(M)-1,1);
-
-for iter = 1:n_skip:max_iter
-
-det_M = det(M(1:iter,1:iter));
-
-        if det_M <= 0
-            not_P = 1;
-            break
-        end
-
-end
-
-    end
-
-end
+    not_P = 1;
 
+end
\ No newline at end of file
diff --git a/Print.m b/Print.m
new file mode 100644
index 0000000..f5bff66
--- /dev/null
+++ b/Print.m
@@ -0,0 +1,12 @@
+%Print.m
+%Messages to print to screen after using the algorithm
+
+if not_P == 0; disp('M matrix is a P-matrix'); 
+elseif not_P == 1; disp('M matrix is not a P-matrix');
+end
+
+if Message_1 == 1
+    x_disp = sprintf('Increase N_guess if you want to include all %d guesses generated by Guesses_master or Guesses_master_2',d0); disp(x_disp)
+elseif Message_2 == 1
+    x_disp = sprintf('More guesses made than are generated by the Guesses_master file used: %d versus %d. The excess guesses are random guesses of 0s and 1s.',N_guess,d0); disp(x_disp)
+end
\ No newline at end of file
diff --git a/Select_solution.m b/Select_solution.m
index 09bc1e2..8823952 100644
--- a/Select_solution.m
+++ b/Select_solution.m
@@ -1,4 +1,4 @@
-%Select a solution (resolving the indeterminacy). Written by Michael Hatcher (m.c.hatcher@soton.ac.uk). 
+%Select a solution (Remark 1, resolving the indeterminacy). Written by Michael Hatcher (m.c.hatcher@soton.ac.uk). 
 %Any errors are my own. Last updated: 12/04/2023.
 
 %Pick a solution (resolving indeterminacy)
diff --git a/Solutions_insert.m b/Solutions_insert.m
new file mode 100644
index 0000000..3df142c
--- /dev/null
+++ b/Solutions_insert.m
@@ -0,0 +1,37 @@
+%Solutions_insert
+
+X_exog = [];
+mstar(mstar==0) = [];
+
+if isempty(mstar) 
+    
+    disp('No solution found. Check T_guess and N_guess are not too small.')
+    
+elseif ~isempty(mstar) 
+
+    X_sol = X_sol(:,:,mstar);
+    X_exog = X_sol(nvar-nx+1:nvar,:,1);
+    X_sol_exc = X_sol_exc(:,:,mstar);   
+
+    solutions =  reshape(permute(X_sol_exc,[1,3,2]),[],size(X_sol_exc,2));
+    ind_solutions = ind_sol(mstar,:);
+
+    %Solution paths
+    x_fin = unique(solutions,'rows','stable');  
+    sol_ind = unique(ind_solutions,'rows','stable');
+
+    %No. of solutions
+    k = length(x_fin(:,1))/(nvar-nx);
+
+    if k==1
+        disp('Unique solution found')
+    else
+        disp('Multiple solutions found')
+    end
+
+    %Final solutions
+    ind_fin = sol_ind; 
+    x_fin = [x_fin; X_exog];  %Put exogenous last
+
+    
+end
diff --git a/Solutions_insert_App_2.m b/Solutions_insert_App_2.m
new file mode 100644
index 0000000..f8ddc13
--- /dev/null
+++ b/Solutions_insert_App_2.m
@@ -0,0 +1,47 @@
+%Solutions_insert_App_2
+
+mstar(mstar==0) = [];
+
+if isempty(mstar) 
+    
+    disp('No solution found. Check T_guess and N_guess are not too small.')
+    
+elseif ~isempty(mstar) 
+
+    X_sol = X_sol(:,:,mstar);
+    X_exog = X_sol(nvar-nx+1:nvar,:,1);
+    X_sol_exc = X_sol_exc(:,:,mstar);   
+
+    %Note: The below is an ad-hoc workaround. Because in the model i(t) = i*(t) in the reference regime (slack), 
+    %some of these paths will removed by the 'unique' function. To prevent this we add small positive number to 
+    %final simulated point of variable 1 (i.e. i(t)).   
+    for ii=1:length(X_sol(1,1,:))
+        rng(1), noise = abs(randn)*1e-18;
+        orig = X_sol(1,T_sim-1,ii); orig_exc = X_sol_exc(1,T_sim-1,ii); 
+        X_sol(1,T_sim-1,ii) = orig + noise/4;
+        X_sol_exc(1,T_sim-1,ii) = orig_exc + noise/4;
+    end  
+
+    solutions =  reshape(permute(X_sol_exc,[1,3,2]),[],size(X_sol_exc,2));
+    ind_solutions = ind_sol(mstar,:);
+
+    %Solution paths
+    x_fin = unique(solutions,'rows','stable');  
+    sol_ind = unique(ind_solutions,'rows','stable');
+
+    %No. of solutions
+    k = length(x_fin(:,1))/(nvar-nx);
+
+    if k==1
+        disp('Unique solution found')
+    else
+        disp('Multiple solutions found')
+    end
+
+    %Final solutions
+    ind_fin = sol_ind; 
+    x_fin = [x_fin; X_exog];  %Put exogenous last
+    
+end
+
+
diff --git a/Solutions_insert_App_2_loop.m b/Solutions_insert_App_2_loop.m
new file mode 100644
index 0000000..dcd9d70
--- /dev/null
+++ b/Solutions_insert_App_2_loop.m
@@ -0,0 +1,65 @@
+%Solutions_insert_App_2_loop
+
+mstar(mstar==0) = [];
+
+if isempty(mstar) 
+    
+        disp('No solution found. Check T_guess and N_guess are not too small.')
+        count_no_sol(j) = 1;
+        count_unique(j) = 0;
+    
+elseif ~isempty(mstar) 
+
+        X_sol = X_sol(:,:,mstar);
+        X_exog = X_sol(nvar-nx+1:nvar,:,1);  %Order exog last in x vector
+        X_sol_exc = X_sol_exc(:,:,mstar); 
+        X_star_sol = X_star_sol(mstar,:); 
+
+    %Note: The below is an ad-hoc workaround. Because in the model i(t) = i*(t) in the reference regime (slack), 
+    %some of these paths will removed by the 'unique' function. To prevent this we add small positive number to 
+    %final simulated point of variable 1 (i.e. i(t)).   
+    for ii=1:length(X_sol(1,1,:))
+        rng(1), noise = abs(randn)*1e-18;
+        orig = X_sol(1,T_sim-1,ii); orig_exc = X_sol_exc(1,T_sim-1,ii); 
+        X_sol(1,T_sim-1,ii) = orig + noise/4;
+        X_sol_exc(1,T_sim-1,ii) = orig_exc + noise/4;
+    end 
+
+    solutions =  reshape(permute(X_sol_exc,[1,3,2]),[],size(X_sol_exc,2));
+    ind_solutions = ind_sol(mstar,:);
+
+    %Solution paths
+    x_fin = unique(solutions,'rows','stable');  
+    sol_ind = unique(ind_solutions,'rows','stable');
+
+    count_no_sol(j) = 0;
+
+    %No. of solutions
+    k = length(x_fin(:,1))/(nvar-nx);
+
+    %Solutions for sim j
+    X_star_fin = unique(X_star_sol,'rows','stable');
+    ind_fin = sol_ind; 
+    x_fin = [x_fin; X_exog];  %Put exogenous last
+
+     if k==1
+        %disp('Unique solution found')
+        count_multi(j) = 0; count_unique(j) = 1;
+        count(j) = sum(ind_fin(1,:)==0);
+     
+     elseif k>=2
+        %disp('Warning: Multiple solutions or no solution')
+        count_multi(j) = 1; 
+            if k==2
+                count(j) = sum(ind_fin(2,:)==0);  %Count periods at bound
+            else
+                disp('Warning: More than two solutions')
+            end
+     
+     end
+
+end
+
+
+
+