diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index adc7f02e4..efac88c3a 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -229,6 +229,9 @@ elseif ismember(options.solve_algo, [2, 12, 4])
     fre = false;
     for i=length(r)-1:-1:1
         blocklength = r(i+1)-r(i);
+        if s(i+1)-s(i) ~= blocklength
+            error('DYNARE_SOLVE (solve_algo=2|4|12): the Dulmage-Mendelsohn decomposition returned a non-square block. This means that the Jacobian is singular. You may want to try another value for solve_algo.')
+        end
         if options.debug
             dprintf('DYNARE_SOLVE (solve_algo=2|4|12): solving block %u of size %u.', i, blocklength);
         end
diff --git a/mex/sources/block_trust_region/mexFunction.f08 b/mex/sources/block_trust_region/mexFunction.f08
index 1092d4233..512865c6e 100644
--- a/mex/sources/block_trust_region/mexFunction.f08
+++ b/mex/sources/block_trust_region/mexFunction.f08
@@ -207,7 +207,9 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
        f_indices => blocks(i)%row_indices
        x_all => x
        if (size(x_indices) /= size(f_indices)) then
-          call mexErrMsgTxt("Non-square block")
+          call mexErrMsgTxt("DYNARE_SOLVE (solve_algo=13|14): the Dulmage-Mendelsohn &
+               &decomposition returned a non-square block. This means that the &
+               &Jacobian is singular. You may want to try another value for solve_algo.")
        end if
        x_block = x(x_indices)
        call trust_region_solve(x_block, matlab_fcn, info, tolx, tolf, maxiter, factor)