diff --git a/matlab/compute_moments_varendo.m b/matlab/compute_moments_varendo.m
index 1009c32ce..294a9b978 100644
--- a/matlab/compute_moments_varendo.m
+++ b/matlab/compute_moments_varendo.m
@@ -145,7 +145,7 @@ if M_.exo_nbr > 1
     end
     skipline();
     if ~all(M_.H==0)
-        if isoctave
+        if isoctave || matlab_ver_less_than('8.1')
             [observable_name_requested_vars, varlist_pos] = intersect_stable(var_list_, options_.varobs);
         else
             [observable_name_requested_vars, varlist_pos] = intersect(var_list_, options_.varobs, 'stable');
diff --git a/matlab/conditional_variance_decomposition.m b/matlab/conditional_variance_decomposition.m
index b618e173d..5f17f7d8f 100644
--- a/matlab/conditional_variance_decomposition.m
+++ b/matlab/conditional_variance_decomposition.m
@@ -89,7 +89,7 @@ end
 % Measurement error
 
 if ~all(StateSpaceModel.measurement_error==0)
-    if isoctave
+    if isoctave || matlab_ver_less_than('8.1')
         [observable_pos,index_subset,index_observables]=intersect_stable(SubsetOfVariables,StateSpaceModel.observable_pos);
     else
         [observable_pos,index_subset,index_observables]=intersect(SubsetOfVariables,StateSpaceModel.observable_pos,'stable');
@@ -105,4 +105,4 @@ if ~all(StateSpaceModel.measurement_error==0)
     ConditionalVarianceDecomposition_ME(:,:,number_of_state_innovations+1)=1-sum(ConditionalVarianceDecomposition_ME(:,:,1:number_of_state_innovations),3);
 else
     ConditionalVarianceDecomposition_ME=[];    
-end
\ No newline at end of file
+end
diff --git a/matlab/conditional_variance_decomposition_ME_mc_analysis.m b/matlab/conditional_variance_decomposition_ME_mc_analysis.m
index c6b2d7e9d..3daf3fc48 100644
--- a/matlab/conditional_variance_decomposition_ME_mc_analysis.m
+++ b/matlab/conditional_variance_decomposition_ME_mc_analysis.m
@@ -63,7 +63,7 @@ if isempty(exogenous_variable_index)
     end
 end
 
-if isoctave
+if isoctave || matlab_ver_less_than('8.1')
     [observable_pos_requested_vars,index_subset,index_observables]=intersect_stable(var_list,options_.varobs);
 else
     [observable_pos_requested_vars,index_subset,index_observables]=intersect(var_list,options_.varobs,'stable');
@@ -136,4 +136,4 @@ oo_.(FirstField).dsge.ConditionalVarianceDecompositionME.HPDsup.(name_1).(name_2
 oo_.(FirstField).dsge.ConditionalVarianceDecompositionME.deciles.(name_1).(name_2)  = p_deciles;
 if options_.estimation.moments_posterior_density.indicator
     oo_.(FirstField).dsge.ConditionalVarianceDecompositionME.density.(name_1).(name_2) = p_density;
-end
\ No newline at end of file
+end
diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m
index bd9744d36..ee62fb3f7 100644
--- a/matlab/disp_moments.m
+++ b/matlab/disp_moments.m
@@ -50,7 +50,7 @@ y = y(ivar,options_.drop+1:end)';
 
 ME_present=0;
 if ~all(M_.H==0)
-    if isoctave
+    if isoctave || matlab_ver_less_than('8.1')
         [observable_pos_requested_vars, index_subset, index_observables] = intersect_stable(ivar, options_.varobs_id);
     else
         [observable_pos_requested_vars, index_subset, index_observables] = intersect(ivar, options_.varobs_id, 'stable');
diff --git a/matlab/disp_th_moments.m b/matlab/disp_th_moments.m
index fd99502d9..7441330cd 100644
--- a/matlab/disp_th_moments.m
+++ b/matlab/disp_th_moments.m
@@ -55,7 +55,7 @@ oo_.var = oo_.gamma_y{1};
 
 ME_present=0;
 if ~all(M_.H==0)
-    if isoctave
+    if isoctave || matlab_ver_less_than('8.1')
         [observable_pos_requested_vars,index_subset,index_observables]=intersect_stable(ivar,options_.varobs_id);
     else
         [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
@@ -105,7 +105,7 @@ if size(stationary_vars, 1) > 0
             lh = cellofchararraymaxlength(M_.endo_names(ivar(stationary_vars)))+2;
             dyntable(options_, title, headers, M_.endo_names(ivar(stationary_vars)), 100*oo_.gamma_y{options_.ar+2}(stationary_vars,:), lh, 8, 2);
             if ME_present
-                if isoctave
+                if isoctave || matlab_ver_less_than('8.1')
                     [stationary_observables, pos_index_subset] = intersect_stable(index_subset, stationary_vars);
                 else
                     [stationary_observables, pos_index_subset] = intersect(index_subset, stationary_vars, 'stable');
diff --git a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m b/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
index 7ff83d28e..21aa0e032 100644
--- a/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
+++ b/matlab/dsge_simulated_theoretical_conditional_variance_decomposition.m
@@ -84,7 +84,7 @@ MaXNumberOfConditionalDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSaved
 
 ME_present=0;
 if ~all(M_.H==0)
-    if isoctave
+    if isoctave || matlab_ver_less_than('8.1')
         [observable_pos_requested_vars,index_subset,index_observables]=intersect_stable(ivar,options_.varobs_id);
     else
         [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
diff --git a/matlab/dsge_simulated_theoretical_variance_decomposition.m b/matlab/dsge_simulated_theoretical_variance_decomposition.m
index 55aa9e0bd..0037ec427 100644
--- a/matlab/dsge_simulated_theoretical_variance_decomposition.m
+++ b/matlab/dsge_simulated_theoretical_variance_decomposition.m
@@ -87,7 +87,7 @@ MaXNumberOfDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPer
 
 ME_present=0;
 if ~all(M_.H==0)
-    if isoctave
+    if isoctave || matlab_ver_less_than('8.1')
         [observable_pos_requested_vars,index_subset,index_observables]=intersect_stable(ivar,options_.varobs_id);
     else
         [observable_pos_requested_vars,index_subset,index_observables]=intersect(ivar,options_.varobs_id,'stable');
diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index 68bf11cd6..a5ec68255 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -93,8 +93,8 @@ if isoctave && octave_ver_less_than('4.4') && ~user_has_octave_forge_package('na
     p{end+1} = '/missing/corrcoef';
 end
 
-% intersect(..., 'stable') doesn't exist in Octave
-if isoctave
+%% intersect(…, 'stable') doesn't exist in Octave and in MATLAB < R2013a
+if isoctave || matlab_ver_less_than('8.1')
     p{end+1} = '/missing/intersect_stable';
 end
 
diff --git a/matlab/missing/intersect_stable/intersect_stable.m b/matlab/missing/intersect_stable/intersect_stable.m
index 71cd3eeb7..be897b6fb 100644
--- a/matlab/missing/intersect_stable/intersect_stable.m
+++ b/matlab/missing/intersect_stable/intersect_stable.m
@@ -1,63 +1,64 @@
-### Crude implementation of insersect(…, 'stable'), which is missing in Octave
+% Crude implementation of intersect(…, 'stable'), which is missing in Octave
 
-### Copyright (C) 2019 Dynare Team
-###
-### This file is part of Dynare.
-###
-### Dynare is free software: you can redistribute it and/or modify
-### it under the terms of the GNU General Public License as published by
-### the Free Software Foundation, either version 3 of the License, or
-### (at your option) any later version.
-###
-### Dynare is distributed in the hope that it will be useful,
-### but WITHOUT ANY WARRANTY; without even the implied warranty of
-### MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-### GNU General Public License for more details.
-###
-### You should have received a copy of the GNU General Public License
-### along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+% Copyright (C) 2019 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 function [c, ia, ib] = intersect_stable(a, b)
 
-    if (nargin != 2)
-        print_usage ();
-    endif
+    if nargin ~= 2
+        error('intersect_stable: needs exactly 2 input arguments');
+    end
 
-    if (isnumeric (a) && isnumeric (b))
+    if isnumeric (a) && isnumeric (b)
         c = [];
-    elseif (iscell (b))
+    elseif iscell (b)
         c = {};
     else
-        c = "";
-    endif
-    ia = ib = [];
+        c = '';
+    end
+    ia = [];
+    ib = [];
 
-    if (isempty (a) || isempty (b))
+    if isempty (a) || isempty (b)
         return
     else
         isrowvec = isrow (a) && isrow (b);
 
         for i = 1:numel(a)
-            if (iscellstr(c))
+            if iscellstr(c)
                 idx = strcmp(a(i), b);
             else
                 idx = a(i) == b;
-            endif
-            if (any(idx) && !ismember(a(i), c))
+            end
+            if any(idx) && ~ismember(a(i), c)
                 c = [c(:); a(i)];
-                if (nargout > 1)
+                if nargout > 1
                     ia = [ia, i];
                     ib = [ib, find(idx)];
-                endif
-            endif
-        endfor
+                end
+            end
+        end
 
-        ## Adjust output orientation for MATLAB compatibility
-        if (isrowvec)
+        %% Adjust output orientation for MATLAB compatibility
+        if isrowvec
             c = c.';
-        endif
-    endif
-endfunction
+        end
+    end
+end
 
 %!test
 %! a = [3 4 1 5];
@@ -80,15 +81,11 @@ endfunction
 %! assert(b(ib), c)
 
 %!test
-%! a = { "defun", "mapcar", "let", "eval-when"};
-%! b = { "setf", "let", "list", "cdr", "defun"};
+%! a = { 'defun', 'mapcar', 'let', 'eval-when'};
+%! b = { 'setf', 'let', 'list', 'cdr', 'defun'};
 %! [c,ia,ib]=intersect_stable(a,b);
-%! assert(c, { "defun", "let" })
+%! assert(c, { 'defun', 'let' })
 %! assert(ia, [1 3])
 %! assert(ib, [5 2])
 %! assert(a(ia), c)
 %! assert(b(ib), c)
-
-## Local variables:
-## mode: Octave
-## End:
diff --git a/matlab/posterior_analysis.m b/matlab/posterior_analysis.m
index f3b76011d..0c8da9def 100644
--- a/matlab/posterior_analysis.m
+++ b/matlab/posterior_analysis.m
@@ -66,7 +66,7 @@ switch type
                                              M_.exo_names,arg2,vartan,arg1,options_.mh_conf_sig,oo_,options_);
     if ~all(M_.H==0)
         if strmatch(arg1,options_.varobs,'exact')
-            if isoctave
+            if isoctave || matlab_ver_less_than('8.1')
                 [observable_name_requested_vars,index_subset,index_observables]=intersect_stable(vartan,options_.varobs);
             else
                 [observable_name_requested_vars,index_subset,index_observables]=intersect(vartan,options_.varobs,'stable');
@@ -97,4 +97,4 @@ switch type
     end        
   otherwise
     disp('Not yet implemented')
-end
\ No newline at end of file
+end