diff --git a/matlab/stochastic_solver/simultxdet.m b/matlab/stochastic_solver/simultxdet.m
index 45512ad55a556d71111efc973f04a8ae56cc1838..07ff89c4acbbc7d7d0119757c958e1f20437c2ef 100644
--- a/matlab/stochastic_solver/simultxdet.m
+++ b/matlab/stochastic_solver/simultxdet.m
@@ -9,14 +9,18 @@ function [y_,int_width,int_width_ME]=simultxdet(y0,ex,ex_det, iorder,var_list,M_
 %    ex_det:    matrix of deterministic exogenous shocks, starting at period 1-M_.maximum_lag
 %    iorder:    order of approximation
 %    var_list:  list of endogenous variables to simulate
-%   int_width_ME:distance between upper bound and
-%                mean forecast when considering measurement error
+%    M_:        Dynare model structure
+%    oo_:       Dynare results structure
+%    options_:  Dynare options structure
+%
 % OUTPUTS:
 %   yf:          mean forecast
 %   int_width:   distance between upper bound and
 %                mean forecast
 %   int_width_ME:distance between upper bound and
 %                mean forecast when considering measurement error
+%   int_width_ME:distance between upper bound and
+%                mean forecast when considering measurement error
 %
 % The forecast horizon is equal to size(ex, 1).
 % The condition size(ex,1)+M_.maximum_lag=size(ex_det,1) must be verified
@@ -123,6 +127,8 @@ elseif iorder == 2
         end
         k1 = k1+1;
     end
+else
+    error('simultxdet.m: order>2 not supported.')
 end
 
 [A,B] = kalman_transition_matrix(dr,nstatic+(1:nspred),1:nc);
@@ -137,7 +143,7 @@ sigma_y = 0;
 
 var_yf=NaN(iter,nvar); %initialize
 for i=1:iter
-    sigma_y1 = ghx1*sigma_y*ghx1'+sigma_u1;
+    sigma_y1 = ghx1*sigma_y*ghx1'+sigma_u1; %only valid at first order, needs to be fixed, see https://git.dynare.org/Dynare/dynare/-/issues/1940
     var_yf(i,:) = diag(sigma_y1)';
     if i == iter
         break