diff --git a/mex/sources/block_trust_region/mexFunction.f08 b/mex/sources/block_trust_region/mexFunction.f08
index 35d9e7d7a..20dad1143 100644
--- a/mex/sources/block_trust_region/mexFunction.f08
+++ b/mex/sources/block_trust_region/mexFunction.f08
@@ -1,4 +1,4 @@
-! Copyright © 2019-2020 Dynare Team
+! Copyright © 2019-2021 Dynare Team
 !
 ! This file is part of Dynare.
 !
@@ -45,43 +45,35 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
 
   if (nrhs < 4 .or. nlhs /= 2) then
      call mexErrMsgTxt("Must have at least 7 inputs and exactly 2 outputs")
-     return
   end if
 
   if (.not. ((mxIsChar(prhs(1)) .and. mxGetM(prhs(1)) == 1) .or. mxIsClass(prhs(1), "function_handle"))) then
      call mexErrMsgTxt("First argument (function) should be a string or a function handle")
-     return
   end if
 
   if (.not. (mxIsDouble(prhs(2)) .and. (mxGetM(prhs(2)) == 1 .or. mxGetN(prhs(2)) == 1))) then
      call mexErrMsgTxt("Second argument (initial guess) should be a real vector")
-     return
   end if
 
   if (.not. (mxIsScalar(prhs(3)) .and. mxIsNumeric(prhs(3)))) then
      call mexErrMsgTxt("Third argument (tolf) should be a numeric scalar")
-     return
   end if
 
   if (.not. (mxIsScalar(prhs(4)) .and. mxIsNumeric(prhs(4)))) then
      call mexErrMsgTxt("Fourth argument (tolx) should be a numeric scalar")
-     return
   end if
 
   if (.not. (mxIsScalar(prhs(5)) .and. mxIsNumeric(prhs(5)))) then
      call mexErrMsgTxt("Fifth argument (maxiter) should be a numeric scalar")
-     return
   end if
 
   if (.not. (mxIsLogicalScalar(prhs(6)))) then
      call mexErrMsgTxt("Sixth argument (debug) should be a logical scalar")
-     return
   end if
 
   if (.not. (mxIsStruct(prhs(7)) .and. &
        (mxGetNumberOfFields(prhs(7)) == 0 .or. mxGetNumberOfFields(prhs(7)) == 4))) then
      call mexErrMsgTxt("Seventh argument should be a struct with either 0 or 4 fields")
-     return
   end if
   specializedunivariateblocks = (mxGetNumberOfFields(prhs(7)) == 4)
 
@@ -101,7 +93,6 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
        tmp = mxGetField(prhs(7), 1_mwIndex, "isloggedlhs")
        if (.not. (c_associated(tmp) .and. mxIsLogical(tmp) .and. mxGetNumberOfElements(tmp) == size(x))) then
           call mexErrMsgTxt("Seventh argument must have a 'isloggedlhs' field of type logical, of same size as second argument")
-          return
        end if
        isloggedlhs => mxGetLogicals(tmp)
 
@@ -109,20 +100,17 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
        if (.not. (c_associated(tmp) .and. mxIsLogical(tmp) .and. mxGetNumberOfElements(tmp) == size(x))) then
           call mexErrMsgTxt("Seventh argument must have a 'isauxdiffloggedrhs' field of type &
                &logical, of same size as second argument")
-          return
        end if
        isauxdiffloggedrhs => mxGetLogicals(tmp)
 
        lhs = mxGetField(prhs(7), 1_mwIndex, "lhs")
        if (.not. (c_associated(lhs) .and. mxIsCell(lhs) .and. mxGetNumberOfElements(lhs) == size(x))) then
           call mexErrMsgTxt("Seventh argument must have a 'lhs' field of type cell, of same size as second argument")
-          return
        end if
 
        endo_names = mxGetField(prhs(7), 1_mwIndex, "endo_names")
        if (.not. (c_associated(endo_names) .and. mxIsCell(endo_names) .and. mxGetNumberOfElements(endo_names) == size(x))) then
           call mexErrMsgTxt("Seventh argument must have a 'endo_names' field of type cell, of same size as second argument")
-          return
        end if
      end block
 
@@ -194,7 +182,6 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
                   end if
                else
                   call mexErrMsgTxt("Algorithm solve_algo=14 cannot be used with this nonlinear problem")
-                  return
                end if
             end if
           end block
@@ -214,7 +201,6 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
        x_all => x
        if (size(x_indices) /= size(f_indices)) then
           call mexErrMsgTxt("Non-square block")
-          return
        end if
        x_block = x(x_indices)
        call trust_region_solve(x_block, matlab_fcn, info, tolx, tolf, maxiter)
diff --git a/mex/sources/disclyap_fast/disclyap_fast.f08 b/mex/sources/disclyap_fast/disclyap_fast.f08
index 7ed5ae4b3..f33ee303f 100644
--- a/mex/sources/disclyap_fast/disclyap_fast.f08
+++ b/mex/sources/disclyap_fast/disclyap_fast.f08
@@ -20,7 +20,7 @@
 ! This is a Fortran translation of a code originally written by Joe Pearlman
 ! and Alejandro Justiniano.
 
-! Copyright © 2020 Dynare Team
+! Copyright © 2020-2021 Dynare Team
 !
 ! This file is part of Dynare.
 !
@@ -61,7 +61,6 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
 
   if (nlhs < 1 .or. nlhs > 2 .or. nrhs < 3 .or. nrhs > 5) then
      call mexErrMsgTxt("disclyap_fast: requires between 3 and 5 input arguments, and 1 or 2 output arguments")
-     return
   end if
 
   n = mxGetM(prhs(1))
@@ -69,19 +68,16 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
       .or. .not. mxIsDouble(prhs(2)) .or. mxIsComplex(prhs(2)) &
       .or. mxGetN(prhs(1)) /= n .or. mxGetM(prhs(2)) /= n .or. mxGetN(prhs(2)) /= n) then
      call mexErrMsgTxt("disclyap_fast: first two arguments should be real matrices of the same dimension")
-     return
   end if
 
   if (.not. (mxIsScalar(prhs(3)) .and. mxIsNumeric(prhs(3)))) then
      call mexErrMsgTxt("disclyap_fast: third argument (tol) should be a numeric scalar")
-     return
   end if
   tol = mxGetScalar(prhs(3))
 
   if (nrhs >= 4) then
      if (.not. (mxIsLogicalScalar(prhs(4)))) then
         call mexErrMsgTxt("disclyap_fast: fourth argument (check_flag) should be a logical scalar")
-        return
      end if
      check_flag = mxGetScalar(prhs(4)) == 1_c_double
   else
@@ -91,7 +87,6 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
   if (nrhs >= 5) then
      if (.not. (mxIsScalar(prhs(5)) .and. mxIsNumeric(prhs(5)))) then
         call mexErrMsgTxt("disclyap_fast: fifth argument (max_iter) should be a numeric scalar")
-        return
      end if
      max_iter = int(mxGetScalar(prhs(5)))
   else
diff --git a/mex/sources/mjdgges/mjdgges.F08 b/mex/sources/mjdgges/mjdgges.F08
index 74de12568..5d67e41fd 100644
--- a/mex/sources/mjdgges/mjdgges.F08
+++ b/mex/sources/mjdgges/mjdgges.F08
@@ -80,7 +80,6 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
 
   if (nrhs < 2 .or. nrhs > 4 .or. nlhs /= 6) then
      call mexErrMsgTxt("MJDGGES: takes 2, 3 or 4 input arguments and exactly 6 output arguments.")
-     return
   end if
 
   n = mxGetM(prhs(1))
@@ -88,14 +87,12 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
       .or. .not. mxIsDouble(prhs(2)) .or. mxIsComplex(prhs(2)) &
       .or. mxGetN(prhs(1)) /= n .or. mxGetM(prhs(2)) /= n .or. mxGetN(prhs(2)) /= n) then
      call mexErrMsgTxt("MJDGGES: first two arguments should be real matrices of the same dimension")
-     return
   end if
 
   ! Set criterium for stable eigenvalues
   if (nrhs >= 3 .and. mxGetM(prhs(3)) > 0) then
      if (.not. (mxIsScalar(prhs(3)) .and. mxIsNumeric(prhs(3)))) then
         call mexErrMsgTxt("MJDGGES: third argument (qz_criterium) should be a numeric scalar")
-        return
      end if
      criterium = mxGetScalar(prhs(3))
   else
@@ -106,7 +103,6 @@ subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
   if (nrhs == 4 .and. mxGetM(prhs(4)) > 0) then
      if (.not. (mxIsScalar(prhs(4)) .and. mxIsNumeric(prhs(4)))) then
         call mexErrMsgTxt("MJDGGES: fourth argument (zhreshold) should be a numeric scalar")
-        return
      end if
      zhreshold = mxGetScalar(prhs(4))
   else