diff --git a/matlab/dr1.m b/matlab/dr1.m
index e9976fd6eefb7f098662717629e7ffdedea7b8ae..d5fbf267bb6b2978d615916e6185861e7ea2d451 100644
--- a/matlab/dr1.m
+++ b/matlab/dr1.m
@@ -555,7 +555,7 @@ zx=[zx; zeros(M_.exo_nbr,np);zeros(M_.exo_det_nbr,np)];
 zu=[zu; eye(M_.exo_nbr);zeros(M_.exo_det_nbr,M_.exo_nbr)];
 [nrzx,nczx] = size(zx);
 
-[err, rhs] = sparse_hessian_times_B_kronecker_C(hessian,zx,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
+[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian,zx,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
 mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
 rhs = -rhs;
 
@@ -615,7 +615,7 @@ hu = dr.ghu(nstatic+1:nstatic+npred,:);
 %kk = kk(1:npred,1:npred);
 %rhs = -hessian*kron(zx,zu)-f1*dr.ghxx(end-nyf+1:end,kk(:))*kron(hx(1:npred,:),hu(1:npred,:));
 
-[err, rhs] = sparse_hessian_times_B_kronecker_C(hessian,zx,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
+[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian,zx,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
 mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
 
 nyf1 = sum(kstate(:,2) == M_.maximum_endo_lag+2);
@@ -624,7 +624,7 @@ hu1 = [hu;zeros(np-npred,M_.exo_nbr)];
 [nrhx,nchx] = size(hx);
 [nrhu1,nchu1] = size(hu1);
 
-[err, abcOut] = A_times_B_kronecker_C(dr.ghxx,hx,hu1,options_.threads.kronecker.A_times_B_kronecker_C);
+[abcOut,err] = A_times_B_kronecker_C(dr.ghxx,hx,hu1,options_.threads.kronecker.A_times_B_kronecker_C);
 mexErrCheck('A_times_B_kronecker_C', err);
 B1 = B*abcOut;
 rhs = -[rhs; zeros(n-M_.endo_nbr,size(rhs,2))]-B1;
@@ -638,10 +638,10 @@ dr.ghxu = A\rhs;
 kk = reshape([1:np*np],np,np);
 kk = kk(1:npred,1:npred);
 
-[err, rhs] = sparse_hessian_times_B_kronecker_C(hessian,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
+[rhs, err] = sparse_hessian_times_B_kronecker_C(hessian,zu,options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
 mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
 
-[err, B1] = A_times_B_kronecker_C(B*dr.ghxx,hu1,options_.threads.kronecker.A_times_B_kronecker_C);
+[B1, err] = A_times_B_kronecker_C(B*dr.ghxx,hu1,options_.threads.kronecker.A_times_B_kronecker_C);
 mexErrCheck('A_times_B_kronecker_C', err);
 rhs = -[rhs; zeros(n-M_.endo_nbr,size(rhs,2))]-B1;
 
@@ -686,7 +686,7 @@ for i=1:M_.maximum_endo_lead
         [junk,k3a,k3] = ...
             find(M_.lead_lag_incidenceordered(M_.maximum_endo_lag+j+1,:));
         nk3a = length(k3a);
-        [err, B1] = sparse_hessian_times_B_kronecker_C(hessian(:,kh(k3,k3)),gu(k3a,:),options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
+        [B1, err] = sparse_hessian_times_B_kronecker_C(hessian(:,kh(k3,k3)),gu(k3a,:),options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
         mexErrCheck('sparse_hessian_times_B_kronecker_C', err);
         RHS = RHS + jacobia_(:,k2)*guu(k2a,:)+B1;
     end
@@ -701,9 +701,9 @@ for i=1:M_.maximum_endo_lead
     kk = find(kstate(:,2) == M_.maximum_endo_lag+i+1);
     gu = dr.ghx*Gu;
     [nrGu,ncGu] = size(Gu);
-    [err, G1] = A_times_B_kronecker_C(dr.ghxx,Gu,options_.threads.kronecker.A_times_B_kronecker_C);
+    [G1, err] = A_times_B_kronecker_C(dr.ghxx,Gu,options_.threads.kronecker.A_times_B_kronecker_C);
     mexErrCheck('A_times_B_kronecker_C', err);
-    [err, G2] = A_times_B_kronecker_C(hxx,Gu,options_.threads.kronecker.A_times_B_kronecker_C);
+    [G2, err] = A_times_B_kronecker_C(hxx,Gu,options_.threads.kronecker.A_times_B_kronecker_C);
     mexErrCheck('A_times_B_kronecker_C', err);
     guu = dr.ghx*Guu+G1;
     Gu = hx*Gu;
diff --git a/matlab/evaluate_planner_objective.m b/matlab/evaluate_planner_objective.m
index b3b8b4e7d1f2938d03f8f2db3bfd88dfd7bf6888..328ba601502feffe4e962c0c220b591932b927b8 100644
--- a/matlab/evaluate_planner_objective.m
+++ b/matlab/evaluate_planner_objective.m
@@ -77,11 +77,11 @@ u = oo.exo_simul(1,:)';
 
 Uyy = full(Uyy);
 
-[err,Uyygygy] = A_times_B_kronecker_C(Uyy,gy,gy,options.threads.kronecker.A_times_B_kronecker_C);
+[Uyygygy, err] = A_times_B_kronecker_C(Uyy,gy,gy,options.threads.kronecker.A_times_B_kronecker_C);
 mexErrCheck('A_times_B_kronecker_C', err);
-[err,Uyygugu] = A_times_B_kronecker_C(Uyy,gu,gu,options.threads.kronecker.A_times_B_kronecker_C);
+[Uyygugu, err] = A_times_B_kronecker_C(Uyy,gu,gu,options.threads.kronecker.A_times_B_kronecker_C);
 mexErrCheck('A_times_B_kronecker_C', err);
-[err,Uyygygu] = A_times_B_kronecker_C(Uyy,gy,gu,options.threads.kronecker.A_times_B_kronecker_C);
+[Uyygygu, err] = A_times_B_kronecker_C(Uyy,gy,gu,options.threads.kronecker.A_times_B_kronecker_C);
 mexErrCheck('A_times_B_kronecker_C', err);
 
 Wbar =U/(1-beta);
@@ -92,9 +92,9 @@ if LQ
 else
     Wyy = (Uy*gyy+Uyygygy+beta*Wy*Gyy)/(eye(npred*npred)-beta*kron(Gy,Gy));
 end
-[err,Wyygugu] = A_times_B_kronecker_C(Wyy,Gu,Gu,options.threads.kronecker.A_times_B_kronecker_C);
+[Wyygugu, err] = A_times_B_kronecker_C(Wyy,Gu,Gu,options.threads.kronecker.A_times_B_kronecker_C);
 mexErrCheck('A_times_B_kronecker_C', err);
-[err,Wyygygu] = A_times_B_kronecker_C(Wyy,Gy,Gu,options.threads.kronecker.A_times_B_kronecker_C);
+[Wyygygu,err] = A_times_B_kronecker_C(Wyy,Gy,Gu,options.threads.kronecker.A_times_B_kronecker_C);
 mexErrCheck('A_times_B_kronecker_C', err);
 if LQ
     Wuu = Uyygugu+beta*Wyygugu;
@@ -114,11 +114,11 @@ end
 yhat = yhat(dr.order_var(nstatic+(1:npred)),1)-dr.ys(dr.order_var(nstatic+(1:npred)));
 u = oo.exo_simul(1,:)';
 
-[err,Wyyyhatyhat] = A_times_B_kronecker_C(Wyy,yhat,yhat,options.threads.kronecker.A_times_B_kronecker_C);
+[Wyyyhatyhat, err] = A_times_B_kronecker_C(Wyy,yhat,yhat,options.threads.kronecker.A_times_B_kronecker_C);
 mexErrCheck('A_times_B_kronecker_C', err);
-[err,Wuuuu] = A_times_B_kronecker_C(Wuu,u,u,options.threads.kronecker.A_times_B_kronecker_C);
+[Wuuuu, err] = A_times_B_kronecker_C(Wuu,u,u,options.threads.kronecker.A_times_B_kronecker_C);
 mexErrCheck('A_times_B_kronecker_C', err);
-[err,Wyuyhatu] = A_times_B_kronecker_C(Wyu,yhat,u,options.threads.kronecker.A_times_B_kronecker_C);
+[Wyuyhatu, err] = A_times_B_kronecker_C(Wyu,yhat,u,options.threads.kronecker.A_times_B_kronecker_C);
 mexErrCheck('A_times_B_kronecker_C', err);
 planner_objective_value(1) = Wbar+Wy*yhat+Wu*u+Wyuyhatu ...
     + 0.5*(Wyyyhatyhat + Wuuuu+Wss);
diff --git a/matlab/kronecker/A_times_B_kronecker_C.m b/matlab/kronecker/A_times_B_kronecker_C.m
index c70c2e86a46232505ab4135c3ca081281d2c9d66..626054ac532d0cb5f1f4ce6f2691fc7101494969 100644
--- a/matlab/kronecker/A_times_B_kronecker_C.m
+++ b/matlab/kronecker/A_times_B_kronecker_C.m
@@ -1,4 +1,39 @@
-function [err, D] = A_times_B_kronecker_C(A,B,C,fake)
+function [D, err] = A_times_B_kronecker_C(A,B,C,fake)
+
+%@info:
+%! @deftypefn {Function File} {@var{D}, @var{err} =} A_timesB_kronecker_C (@var{A},@var{B},@var{C},@var{fake})
+%! @anchor{kronecker/A_times_B_kronecker_C}
+%! @sp 1
+%! Computes A*kron(B,C).
+%! @sp 2
+%! @strong{Inputs}
+%! @sp 1
+%! @table @ @var
+%! @item A
+%! mA*nA matrix of doubles.
+%! @item B
+%! mB*nB matrix of doubles.
+%! @item C
+%! mC*nC matrix of doubles.
+%! @item fake
+%! Anything you want, just a fake parameter (because the mex version admits a last argument specifying the number of threads to be used in parallel mode).
+%! @end table
+%! @sp 2
+%! @strong{Outputs}
+%! @sp 1
+%! @table @ @var
+%! @item c
+%! Integer scalar, the number of years, quarters, months or weeks between @var{a} and @var{B}.
+%! @end table
+%! @sp 2
+%! @strong{This function is called by:}
+%! @sp 2
+%! @strong{This function calls:}
+%! @ref{@@dynDates/eq},@ref{@@dynDates/lt}
+%!
+%! @end deftypefn
+%@eod:
+
 %function [err, D] = A_times_B_kronecker_C(A,B,C)
 % Computes A * kron(B,C). 
 %
diff --git a/matlab/kronecker/sparse_hessian_times_B_kronecker_C.m b/matlab/kronecker/sparse_hessian_times_B_kronecker_C.m
index 028ea92441405d8f650aa6bcac7c1cd1a590fc43..555834510c484f6bea3cbb8f25b9c1185dc26198 100644
--- a/matlab/kronecker/sparse_hessian_times_B_kronecker_C.m
+++ b/matlab/kronecker/sparse_hessian_times_B_kronecker_C.m
@@ -1,4 +1,4 @@
-function [err, D] = sparse_hessian_times_B_kronecker_C(varargin)
+function [D, err] = sparse_hessian_times_B_kronecker_C(varargin)
 %function [err, D] = sparse_hessian_times_B_kronecker_C(A,B,C, fake)
 % Computes A * kron(B,C) where A is a sparse matrix.
 %
@@ -43,9 +43,9 @@ end
 
 switch nargin
   case 4
-    [fake,D] = A_times_B_kronecker_C(A,B,C,fake);
+    [D, fake] = A_times_B_kronecker_C(A,B,C,fake);
   case 3
-    [fake,D] = A_times_B_kronecker_C(A,B,C);
+    [D, fake] = A_times_B_kronecker_C(A,B,C);
   otherwise
     error('Two or Three input arguments required!')
 end
diff --git a/matlab/simult_.m b/matlab/simult_.m
index 199b354327b19b780356ff8e4e5705618d61f8f4..46823f87e5d4a5dedb496184aa2b6a3ba8a0014b 100644
--- a/matlab/simult_.m
+++ b/matlab/simult_.m
@@ -105,11 +105,11 @@ else
                 yhat1 = y__(order_var(k2))-dr.ys(order_var(k2));
                 yhat2 = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
                 epsilon = ex_(i-1,:)';
-                [err, abcOut1] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1,options_.threads.kronecker.A_times_B_kronecker_C);
+                [abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat1,options_.threads.kronecker.A_times_B_kronecker_C);
                 mexErrCheck('A_times_B_kronecker_C', err);
-                [err, abcOut2] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
+                [abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
                 mexErrCheck('A_times_B_kronecker_C', err);
-                [err, abcOut3] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
+                [abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat1,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
                 mexErrCheck('A_times_B_kronecker_C', err);
                 y_(order_var,i) = constant + dr.ghx*yhat2 + dr.ghu*epsilon ...
                     + abcOut1 + abcOut2 + abcOut3;
@@ -119,11 +119,11 @@ else
             for i = 2:iter+M_.maximum_lag
                 yhat = y_(order_var(k2),i-1)-dr.ys(order_var(k2));
                 epsilon = ex_(i-1,:)';
-                [err, abcOut1] = A_times_B_kronecker_C(.5*dr.ghxx,yhat,options_.threads.kronecker.A_times_B_kronecker_C);
+                [abcOut1, err] = A_times_B_kronecker_C(.5*dr.ghxx,yhat,options_.threads.kronecker.A_times_B_kronecker_C);
                 mexErrCheck('A_times_B_kronecker_C', err);
-                [err, abcOut2] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
+                [abcOut2, err] = A_times_B_kronecker_C(.5*dr.ghuu,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
                 mexErrCheck('A_times_B_kronecker_C', err);
-                [err, abcOut3] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
+                [abcOut3, err] = A_times_B_kronecker_C(dr.ghxu,yhat,epsilon,options_.threads.kronecker.A_times_B_kronecker_C);
                 mexErrCheck('A_times_B_kronecker_C', err);
                 y_(dr.order_var,i) = constant + dr.ghx*yhat + dr.ghu*epsilon ...
                     + abcOut1 + abcOut2 + abcOut3;
diff --git a/mex/sources/kronecker/A_times_B_kronecker_C.cc b/mex/sources/kronecker/A_times_B_kronecker_C.cc
index e30d9d1bf36906ab1c2beb63d0f0f327552ac9b1..901d68b95f768a7135330ca4ba79f8691eef79c0 100644
--- a/mex/sources/kronecker/A_times_B_kronecker_C.cc
+++ b/mex/sources/kronecker/A_times_B_kronecker_C.cc
@@ -126,8 +126,8 @@ void
 mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
   // Check input and output:
-  if (nrhs > 4 || nrhs < 3 || nlhs != 2)
-    DYN_MEX_FUNC_ERR_MSG_TXT("A_times_B_kronecker_C takes 3 or 4 input arguments and provides exactly 2 output arguments.");
+  if (nrhs > 4 || nrhs < 3)
+    DYN_MEX_FUNC_ERR_MSG_TXT("A_times_B_kronecker_C takes 3 or 4 input arguments and provides 2 output arguments.");
 
   // Get & Check dimensions (columns and rows):
   mwSize mA, nA, mB, nB, mC, nC;
@@ -164,13 +164,13 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   double *D;
   if (nrhs == 4)
     {
-      plhs[1] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
+      plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
     }
   else
     {
-      plhs[1] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
+      plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
     }
-  D = mxGetPr(plhs[1]);
+  D = mxGetPr(plhs[0]);
   // Computational part:
   if (nrhs == 3)
     {
@@ -180,5 +180,5 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     {
       full_A_times_kronecker_B_C(A, B, C, &D[0], mA, nA, mB, nB, mC, nC, numthreads);
     }
-  plhs[0] = mxCreateDoubleScalar(0);
+  plhs[1] = mxCreateDoubleScalar(0);
 }
diff --git a/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc b/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc
index 1ba6b8c7f0612da11e9cbd69fc5ef06e25f590e6..d4f50cc57d6096d4fb42bf4045cdc57cc5a5a84e 100644
--- a/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc
+++ b/mex/sources/kronecker/sparse_hessian_times_B_kronecker_C.cc
@@ -143,8 +143,8 @@ void
 mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 {
   // Check input and output:
-  if ((nrhs > 4) || (nrhs < 3) || nlhs != 2)
-    DYN_MEX_FUNC_ERR_MSG_TXT("sparse_hessian_times_B_kronecker_C takes 3 or 4 input arguments and provides exactly 2 output argument.");
+  if ((nrhs > 4) || (nrhs < 3) )
+    DYN_MEX_FUNC_ERR_MSG_TXT("sparse_hessian_times_B_kronecker_C takes 3 or 4 input arguments and provides 2 output arguments.");
 
   if (!mxIsSparse(prhs[0]))
     DYN_MEX_FUNC_ERR_MSG_TXT("sparse_hessian_times_B_kronecker_C: First input must be a sparse (dynare) hessian matrix.");
@@ -185,13 +185,13 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
   double *D;
   if (nrhs == 4)
     {
-      plhs[1] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
+      plhs[0] = mxCreateDoubleMatrix(mA, nB*nC, mxREAL);
     }
   else
     {
-      plhs[1] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
+      plhs[0] = mxCreateDoubleMatrix(mA, nB*nB, mxREAL);
     }
-  D = mxGetPr(plhs[1]);
+  D = mxGetPr(plhs[0]);
   // Computational part:
   if (nrhs == 3)
     {
@@ -201,5 +201,5 @@ mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
     {
       sparse_hessian_times_B_kronecker_C(isparseA, jsparseA, vsparseA, B, C, D, mA, nA, mB, nB, mC, nC, numthreads);
     }
-  plhs[0] = mxCreateDoubleScalar(0);
+  plhs[1] = mxCreateDoubleScalar(0);
 }