diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index f68d14cce3d26c6dfb7597c478417591c97180b5..9e500f59b67b1f7954e0063273805ff346b65f7f 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -370,7 +370,7 @@ elseif ismember(options.solve_algo, [13, 14])
         auxstruct.isloggedlhs = isloggedlhs;
         auxstruct.isauxdiffloggedrhs = isauxdiffloggedrhs;
     end
-    [x, errorflag] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, auxstruct, arguments{:});
+    [x, errorflag, exitflag] = block_trust_region(f, x, tolf, options.solve_tolx, maxit, options.trust_region_initial_step_bound_factor, options.debug, auxstruct, arguments{:});
     [fvec, fjac] = feval(f, x, arguments{:});
 else
     error('DYNARE_SOLVE: option solve_algo must be one of [0,1,2,3,4,9,10,11,12,13,14]')
diff --git a/mex/sources/block_trust_region/mexFunction.f08 b/mex/sources/block_trust_region/mexFunction.f08
index 62398bc300ddf0db20435dafdc0663a69470432d..00e9b73280d90c30ae45bc53809b4f4a7510aa40 100644
--- a/mex/sources/block_trust_region/mexFunction.f08
+++ b/mex/sources/block_trust_region/mexFunction.f08
@@ -43,8 +43,8 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
   logical :: fre ! True if the last block has been solved (i.e. not evaluated), so that residuals must be updated
   integer, dimension(:), allocatable :: evaled_cols ! If fre=.false., lists the columns that have been evaluated so far without updating the residuals
 
-  if (nrhs < 4 .or. nlhs /= 2) then
-     call mexErrMsgTxt("Must have at least 7 inputs and exactly 2 outputs")
+  if (nrhs < 4 .or. nlhs /= 3) then
+     call mexErrMsgTxt("Must have at least 7 inputs and exactly 3 outputs")
   end if
 
   if (.not. ((mxIsChar(prhs(1)) .and. mxGetM(prhs(1)) == 1) .or. mxIsClass(prhs(1), "function_handle"))) then
@@ -231,14 +231,18 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
           call mexPrintf_trim_newline("DYNARE_SOLVE (solve_algo=13|14): residuals still too large, solving for the whole model")
      call trust_region_solve(x, matlab_fcn, info, tolx, tolf, maxiter, factor)
   else
-     info = 1
+     if (size(blocks).gt.1) then
+        ! Note that the value of info may be different across blocks
+        info = 1
+     end if
   end if
 
   plhs(1) = mxCreateDoubleMatrix(int(size(x, 1), mwSize), 1_mwSize, mxREAL)
   mxGetPr(plhs(1)) = x
-  if (info == 1) then
+  if ((info == 1) .or. (info == -1)) then
      plhs(2) = mxCreateDoubleScalar(0._c_double)
   else
      plhs(2) = mxCreateDoubleScalar(1._c_double)
   end if
+  plhs(3) = mxCreateDoubleScalar(dble(info))
 end subroutine mexFunction
diff --git a/tests/nonlinearsolvers.m b/tests/nonlinearsolvers.m
index 5231dcd7b02465d8cbd8c2a29ea9b8f5898a9d90..8c4cfc1e6829985e56e522413c57747e89080c0e 100644
--- a/tests/nonlinearsolvers.m
+++ b/tests/nonlinearsolvers.m
@@ -74,7 +74,7 @@ for i=1:length(objfun)
         x = objfun{i}();
     end
     try
-        [x, errorflag] = block_trust_region(objfun{i}, x, tolf, tolx, maxit, factor, false, auxstruct);
+        [x, errorflag, exitflag] = block_trust_region(objfun{i}, x, tolf, tolx, maxit, factor, false, auxstruct);
         if isequal(func2str(objfun{i}), 'powell2')
             if ~errorflag
                 testFailed = testFailed+1;