diff --git a/mex/build/dynare_minreal.am b/mex/build/dynare_minreal.am
new file mode 100644
index 0000000000000000000000000000000000000000..a0514ecfc6d9ea5144f36a28c9b914372e45a1bd
--- /dev/null
+++ b/mex/build/dynare_minreal.am
@@ -0,0 +1,14 @@
+mex_PROGRAMS = dynare_minreal
+
+dynare_minreal_LDADD = $(LIBADD_SLICOT)
+dynare_minreal_LDFLAGS = $(AM_LDFLAGS) $(LDFLAGS_SLICOT)
+
+nodist_dynare_minreal_SOURCES = dynare_minreal.f08 matlab_mex.F08 slicot.F08
+
+BUILT_SOURCES = $(nodist_dynare_minreal_SOURCES)
+CLEANFILES = $(nodist_dynare_minreal_SOURCES)
+
+dynare_minreal.o : matlab_mex.mod slicot.mod
+
+%.f08: $(top_srcdir)/../../sources/dynare_minreal/%.f08
+	$(LN_S) -f $< $@
diff --git a/mex/build/matlab/Makefile.am b/mex/build/matlab/Makefile.am
index e84baee6ed388fe01510b7a2d7a63d0e045b98b4..f1a7654b8e32af9e608befa0a5f0e0b591a403e4 100644
--- a/mex/build/matlab/Makefile.am
+++ b/mex/build/matlab/Makefile.am
@@ -11,8 +11,8 @@ if ENABLE_MEX_MS_SBVAR
 SUBDIRS += ms_sbvar
 endif
 
-if ENABLE_MEX_KALMAN_STEADY_STATE
-SUBDIRS += kalman_steady_state
+if ENABLE_MEX_KALMAN_MINREAL
+SUBDIRS += kalman_steady_state dynare_minreal
 endif
 
 if HAVE_M2HTML
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index c9460676a4172feab47f6c316d34f1c8f1b22cb1..7df7907c2193712fb4b3bf02b0b454ccaea95a81 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -1,6 +1,6 @@
 dnl Process this file with autoconf to produce a configure script.
 
-dnl Copyright © 2009-2020 Dynare Team
+dnl Copyright © 2009-2021 Dynare Team
 dnl
 dnl This file is part of Dynare.
 dnl
@@ -72,8 +72,8 @@ AM_CONDITIONAL([ENABLE_MEX_MS_SBVAR], [test "$enable_mex_ms_sbvar" = yes])
 AC_ARG_ENABLE([mex-dynare++], AS_HELP_STRING([--disable-mex-dynare++], [disable compilation of MEX based on Dynare++]), [], [enable_mex_dynareplusplus=yes])
 AM_CONDITIONAL([ENABLE_MEX_DYNAREPLUSPLUS], [test "$enable_mex_dynareplusplus" = yes])
 
-AC_ARG_ENABLE([mex-kalman-steady-state], AS_HELP_STRING([--disable-mex-kalman-steady-state], [disable compilation of the kalman_steady_state MEX]), [], [enable_mex_kalman_steady_state=yes])
-AM_CONDITIONAL([ENABLE_MEX_KALMAN_STEADY_STATE], [test "$enable_mex_kalman_steady_state" = yes])
+AC_ARG_ENABLE([mex-kalman-minreal], AS_HELP_STRING([--disable-mex-kalman-minreal], [disable compilation of the kalman_steady_state and dynare_minreal MEX]), [], [enable_mex_kalman_minreal=yes])
+AM_CONDITIONAL([ENABLE_MEX_KALMAN_MINREAL], [test "$enable_mex_kalman_minreal" = yes])
 
 # Check for GSL, needed by MS-SBVAR MEX
 if test "$enable_mex_ms_sbvar" = yes; then
@@ -87,10 +87,10 @@ if test "$enable_mex_dynareplusplus" = yes; then
   test "$has_matio" != yes && AC_MSG_ERROR([libmatio cannot be found. If you want to skip the compilation of MEX files based Dynare++, pass the --disable-mex-dynare++ flag.])
 fi
 
-# Check for libslicot, needed by kalman_steady_state
-if test "$enable_mex_kalman_steady_state" = yes; then
+# Check for libslicot, needed by kalman_steady_state and dynare_minreal
+if test "$enable_mex_kalman_minreal" = yes; then
   AX_SLICOT([matlab])
-  test "$has_slicot" != yes && AC_MSG_ERROR([slicot cannot be found. If you want to skip the compilation of the kalman_steady_state MEX, pass the --disable-mex-kalman-steady-state flag.])
+  test "$has_slicot" != yes && AC_MSG_ERROR([slicot cannot be found. If you want to skip the compilation of the kalman_steady_state and dynare_minreal MEX, pass the --disable-mex-kalman-minreal flag.])
 fi
 
 # On Windows, we want static linking of the external libraries
@@ -120,10 +120,10 @@ else
    BUILD_GENSYLV_KORDER_DYNSIMUL_MEX_MATLAB="no"
 fi
 
-if test "$enable_mex_kalman_steady_state" = yes; then
-   BUILD_KALMAN_STEADY_STATE_MATLAB="yes"
+if test "$enable_mex_kalman_minreal" = yes; then
+   BUILD_KALMAN_MINREAL_MATLAB="yes"
 else
-   BUILD_KALMAN_STEADY_STATE_MATLAB="no"
+   BUILD_KALMAN_MINREAL_MATLAB="no"
 fi
 
 if test "$enable_mex_ms_sbvar" = yes; then
@@ -140,7 +140,7 @@ Binaries (with "make"):
  MEX files for MATLAB (except those listed below):                   yes
  Gensylv, k-order and dynare_simul MEX files for MATLAB:             $BUILD_GENSYLV_KORDER_DYNSIMUL_MEX_MATLAB
  MS-SBVAR MEX files for MATLAB:                                      $BUILD_MS_SBVAR_MEX_MATLAB
- Kalman Steady State MEX file for MATLAB:                            $BUILD_KALMAN_STEADY_STATE_MATLAB
+ Kalman Steady State and dynare_minreal MEX file for MATLAB:         $BUILD_KALMAN_MINREAL_MATLAB
  M2HTML documentation:                                               $BUILD_M2HTML
 
 ])
@@ -161,6 +161,7 @@ AC_CONFIG_FILES([Makefile
                  perfect_foresight_problem/Makefile
                  num_procs/Makefile
                  block_trust_region/Makefile
-                 disclyap_fast/Makefile])
+                 disclyap_fast/Makefile
+                 dynare_minreal/Makefile])
 
 AC_OUTPUT
diff --git a/mex/build/matlab/dynare_minreal/Makefile.am b/mex/build/matlab/dynare_minreal/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..c212c7afaed4d00299b6ac6c786c683c9da7f0ab
--- /dev/null
+++ b/mex/build/matlab/dynare_minreal/Makefile.am
@@ -0,0 +1,2 @@
+include ../mex.am
+include ../../dynare_minreal.am
diff --git a/mex/build/matlab/mex.am b/mex/build/matlab/mex.am
index 5ba01b72306f8c083666f8f115302fe7d4c23249..4d0c52753dab1b0a1307ff93d8ead6d52ffedbf6 100644
--- a/mex/build/matlab/mex.am
+++ b/mex/build/matlab/mex.am
@@ -29,7 +29,7 @@ clean-local:
 		cd $(top_srcdir)/../../matlab && rm -f $(PROGRAMS); \
 	fi
 
-# Rules for the Fortran 2008 interface to MEX and BLAS/LAPACK functions
+# Rules for the Fortran 2008 interface to MEX and BLAS/LAPACK/SLICOT functions
 matlab_mat.mod: matlab_mex.o
 matlab_mex.mod: matlab_mex.o
 
@@ -41,3 +41,8 @@ lapack.mod: blas_lapack.o
 
 blas_lapack.F08: $(top_srcdir)/../../sources/blas_lapack.F08
 	$(LN_S) -f $< $@
+
+slicot.mod: slicot.o
+
+slicot.F08: $(top_srcdir)/../../sources/slicot.F08
+	$(LN_S) -f $< $@
diff --git a/mex/build/octave/Makefile.am b/mex/build/octave/Makefile.am
index 6726be71e928d36dc9cac5cb8deebed59ebf1b60..56c60643add8bbb0f75a22e50503b2f1300cc0b7 100644
--- a/mex/build/octave/Makefile.am
+++ b/mex/build/octave/Makefile.am
@@ -11,8 +11,8 @@ if ENABLE_MEX_MS_SBVAR
 SUBDIRS += ms_sbvar
 endif
 
-if ENABLE_MEX_KALMAN_STEADY_STATE
-SUBDIRS += kalman_steady_state
+if ENABLE_MEX_KALMAN_MINREAL
+SUBDIRS += kalman_steady_state dynare_minreal
 endif
 
 install-exec-local:
diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac
index 6f7f9391549c2f5a734ad47b57fcf70750778469..6592accac9c62b7a8e66c5b4c3a373f86e901bf4 100644
--- a/mex/build/octave/configure.ac
+++ b/mex/build/octave/configure.ac
@@ -1,6 +1,6 @@
 dnl Process this file with autoconf to produce a configure script.
 
-dnl Copyright © 2009-2020 Dynare Team
+dnl Copyright © 2009-2021 Dynare Team
 dnl
 dnl This file is part of Dynare.
 dnl
@@ -70,8 +70,8 @@ AM_CONDITIONAL([ENABLE_MEX_MS_SBVAR], [test "$enable_mex_ms_sbvar" = yes])
 AC_ARG_ENABLE([mex-dynare++], AS_HELP_STRING([--disable-mex-dynare++], [disable compilation of MEX based on Dynare++]), [], [enable_mex_dynareplusplus=yes])
 AM_CONDITIONAL([ENABLE_MEX_DYNAREPLUSPLUS], [test "$enable_mex_dynareplusplus" = yes])
 
-AC_ARG_ENABLE([mex-kalman-steady-state], AS_HELP_STRING([--disable-mex-kalman-steady-state], [disable compilation of the kalman_steady_state MEX]), [], [enable_mex_kalman_steady_state=yes])
-AM_CONDITIONAL([ENABLE_MEX_KALMAN_STEADY_STATE], [test "$enable_mex_kalman_steady_state" = yes])
+AC_ARG_ENABLE([mex-kalman-minreal], AS_HELP_STRING([--disable-mex-kalman-minreal], [disable compilation of the kalman_steady_state and dynare_minreal MEX]), [], [enable_mex_kalman_minreal=yes])
+AM_CONDITIONAL([ENABLE_MEX_KALMAN_MINREAL], [test "$enable_mex_kalman_minreal" = yes])
 
 # Check for GSL, needed by MS-SBVAR
 if test "$enable_mex_ms_sbvar" = yes; then
@@ -85,10 +85,10 @@ if test "$enable_mex_dynareplusplus" = yes -o "$enable_mex_ms_sbvar" = yes; then
   test "$has_matio" != yes && AC_MSG_ERROR([libmatio cannot be found. If you want to skip the compilation of MS-SBVAR MEX and MEX files based Dynare++, pass the --disable-mex-dynare++ and --disable-mex-ms-sbvar flags.])
 fi
 
-# Check for libslicot, needed by kalman_steady_state
-if test "$enable_mex_kalman_steady_state" = yes; then
+# Check for libslicot, needed by kalman_steady_state and dynare_minreal
+if test "$enable_mex_kalman_minreal" = yes; then
   AX_SLICOT([octave])
-  test "$has_slicot" != yes && AC_MSG_ERROR([slicot cannot be found. If you want to skip the compilation of the kalman_steady_state MEX, pass the --disable-mex-kalman-steady-state flag.])
+  test "$has_slicot" != yes && AC_MSG_ERROR([slicot cannot be found. If you want to skip the compilation of the kalman_steady_state and dynare_minreal MEX, pass the --disable-mex-kalman-minreal flag.])
 fi
 
 # Check for UMFPACK, needed by bytecode
@@ -111,10 +111,10 @@ else
    BUILD_GENSYLV_KORDER_DYNSIMUL_MEX_OCTAVE="no"
 fi
 
-if test "$enable_mex_kalman_steady_state" = yes; then
-   BUILD_KALMAN_STEADY_STATE_OCTAVE="yes"
+if test "$enable_mex_kalman_minreal" = yes; then
+   BUILD_KALMAN_MINREAL_OCTAVE="yes"
 else
-   BUILD_KALMAN_STEADY_STATE_OCTAVE="no"
+   BUILD_KALMAN_MINREAL_OCTAVE="no"
 fi
 
 if test "$enable_mex_ms_sbvar" = yes; then
@@ -131,7 +131,7 @@ Binaries (with "make"):
  MEX files for Octave (except those listed below):                   yes
  Gensylv, k-order and dynare_simul MEX for Octave:                   $BUILD_GENSYLV_KORDER_DYNSIMUL_MEX_OCTAVE
  MS-SBVAR MEX files for Octave:                                      $BUILD_MS_SBVAR_MEX_OCTAVE
- Kalman Steady State MEX file for Octave:                            $BUILD_KALMAN_STEADY_STATE_OCTAVE
+ Kalman Steady State and dynare_minreal MEX file for Octave:         $BUILD_KALMAN_MINREAL_OCTAVE
 
 ])
 
@@ -151,6 +151,7 @@ AC_CONFIG_FILES([Makefile
                  perfect_foresight_problem/Makefile
                  num_procs/Makefile
                  block_trust_region/Makefile
-                 disclyap_fast/Makefile])
+                 disclyap_fast/Makefile
+                 dynare_minreal/Makefile])
 
 AC_OUTPUT
diff --git a/mex/build/octave/dynare_minreal/Makefile.am b/mex/build/octave/dynare_minreal/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..c4d28f1706e3f8f29c27dc349cc0d1bfe4b64cf2
--- /dev/null
+++ b/mex/build/octave/dynare_minreal/Makefile.am
@@ -0,0 +1,3 @@
+EXEEXT = .mex
+include ../mex.am
+include ../../dynare_minreal.am
diff --git a/mex/build/octave/mex.am b/mex/build/octave/mex.am
index 77eccb3e7df62bc644d46f4e619ef933dc10b922..41a4f13c246ef02da2e52ee5d72960df8ab919c7 100644
--- a/mex/build/octave/mex.am
+++ b/mex/build/octave/mex.am
@@ -37,7 +37,7 @@ clean-local:
 		cd $(top_srcdir)/../../octave && rm -f $(PROGRAMS); \
 	fi
 
-# Rules for the Fortran 2008 interface to MEX and BLAS/LAPACK functions
+# Rules for the Fortran 2008 interface to MEX and BLAS/LAPACK/SLICOT functions
 matlab_mat.mod: matlab_mex.o
 matlab_mex.mod: matlab_mex.o
 
@@ -49,3 +49,8 @@ lapack.mod: blas_lapack.o
 
 blas_lapack.F08: $(top_srcdir)/../../sources/blas_lapack.F08
 	$(LN_S) -f $< $@
+
+slicot.mod: slicot.o
+
+slicot.F08: $(top_srcdir)/../../sources/slicot.F08
+	$(LN_S) -f $< $@
diff --git a/mex/sources/dynare_minreal/dynare_minreal.f08 b/mex/sources/dynare_minreal/dynare_minreal.f08
new file mode 100644
index 0000000000000000000000000000000000000000..504f9f931b6ecf51519633cf2f4fddd400391d59
--- /dev/null
+++ b/mex/sources/dynare_minreal/dynare_minreal.f08
@@ -0,0 +1,164 @@
+! Equivalent to minreal function from MATLAB’s control toolbox
+! Essentially a wrapper around SLICOT’s TB01PD.
+!
+! Takes as first argument a structure with 4 fields:
+! — A: n×n real matrix
+! − B: n×m real matrix
+! − C: p×n real matrix
+! − D: p×m real matrix
+!
+! Returns a similar structure with 4 fields:
+! — A: nr×nr real matrix
+! − B: nr×m real matrix
+! − C: p×nr real matrix
+! − D: p×m real matrix
+! where nr < n is the size of the minimal state-space system.
+!
+! An optional second argument sets the tolerance (if omitted, a default value
+! is used).
+
+! Copyright © 2021 Dynare Team
+!
+! This file is part of Dynare.
+!
+! Dynare is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! Dynare is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+subroutine mexFunction(nlhs, plhs, nrhs, prhs) bind(c, name='mexFunction')
+  use iso_fortran_env
+  use matlab_mex
+  use slicot
+  implicit none
+
+  type(c_ptr), dimension(*), intent(in), target :: prhs
+  type(c_ptr), dimension(*), intent(out) :: plhs
+  integer(c_int), intent(in), value :: nlhs, nrhs
+
+  real(real64) :: tol
+
+  integer(c_size_t) :: n, m, p
+  integer(slint) :: n_sl, m_sl, p_sl, lda, ldb, ldc, nr, ldwork, info
+  type(c_ptr) :: A_mx, B_mx, C_mx, D_mx, Ar_mx, Br_mx, Cr_mx, Dr_mx
+  real(real64), dimension(:, :), pointer :: A, B, C, Ar, Br, Cr
+  real(real64), dimension(:, :), allocatable :: Atmp, Btmp, Ctmp
+  integer(slint), dimension(:), allocatable :: iwork
+  real(real64), dimension(:), allocatable :: dwork
+  integer(c_int) :: ignore
+
+
+  !! Process and check arguments
+
+  if (nlhs /= 1 .or. nrhs > 2 .or. nrhs < 1) then
+     call mexErrMsgTxt("dynare_minreal: requires 1 or 2 input arguments, and 1 output arguments")
+     return
+  end if
+
+  if (.not. (mxIsStruct(prhs(1)))) then
+     call mexErrMsgTxt("dynare_minreal: first argument (sys) should be a struct")
+     return
+  end if
+
+  A_mx = mxGetField(prhs(1), 1_mwIndex, "A")
+  if (.not. (c_associated(A_mx) .and. mxIsDouble(A_mx) .and. mxGetM(A_mx) == mxGetN(A_mx))) then
+     call mexErrMsgTxt("dynare_minreal: first argument (sys) should have a field A containing a square real matrix")
+     return
+  end if
+  n = mxGetM(A_mx)
+  A(1:n, 1:n) => mxGetPr(A_mx)
+
+  B_mx = mxGetField(prhs(1), 1_mwIndex, "B")
+  if (.not. (c_associated(B_mx) .and. mxIsDouble(B_mx) .and. mxGetM(B_mx) == n)) then
+     call mexErrMsgTxt("dynare_minreal: first argument (sys) should have a field B containing a real &
+          &matrix with as many lines as sys.A")
+     return
+  end if
+  m = mxGetN(B_mx)
+  B(1:n, 1:m) => mxGetPr(B_mx)
+
+  C_mx = mxGetField(prhs(1), 1_mwIndex, "C")
+  if (.not. (c_associated(C_mx) .and. mxIsDouble(C_mx) .and. mxGetN(C_mx) == n)) then
+     call mexErrMsgTxt("dynare_minreal: first argument (sys) should have a field C containing a real &
+          &matrix with as many columns as sys.A")
+     return
+  end if
+  p = mxGetM(C_mx)
+  C(1:p, 1:n) => mxGetPr(C_mx)
+
+  D_mx = mxGetField(prhs(1), 1_mwIndex, "D")
+  if (.not. (c_associated(D_mx) .and. mxIsDouble(D_mx) .and. mxGetM(D_mx) == p .and. mxGetN(D_mx) == m)) then
+     call mexErrMsgTxt("dynare_minreal: first argument (sys) should have a field D containing a real &
+          &matrix with as many rows as sys.C and as many columns as sys.B")
+     return
+  end if
+
+  if (nrhs == 2) then
+      if (.not. (mxIsScalar(prhs(2)) .and. mxIsNumeric(prhs(2)))) then
+         call mexErrMsgTxt("dynare_minreal: second argument (tol) should be a numeric scalar")
+         return
+      end if
+      tol = mxGetScalar(prhs(2))
+   else
+      tol = 0._real64 ! Let SLICOT determine the tolerance
+   end if
+
+
+   !! Construct arguments to be passed to the SLICOT function
+
+   n_sl = int(n, slint)
+   m_sl = int(m, slint)
+   p_sl = int(p, slint)
+   lda = max(1_slint, n_sl)
+   ldb = max(1_slint, n_sl)
+   if (n_sl == 0_slint) then
+      ldc = 1_slint
+   else
+      ldc = max(1_slint, max(m_sl, p_sl))
+   end if
+   ldwork = max(1_slint, n_sl + max(n_sl, max(3_slint*m_sl, 3_slint*p_sl)))
+
+   allocate(Atmp(n_sl, n_sl), Btmp(ldb, max(m_sl, p_sl)), Ctmp(ldc, n_sl))
+   Atmp(1:n_sl, 1:n_sl) = A
+   Btmp(1:n_sl, 1:m_sl) = B
+   Ctmp(1:p_sl, 1:n_sl) = C
+
+   allocate(iwork(n_sl+max(m_sl,p_sl)), dwork(ldwork))
+
+   call tb01pd("M", "N", n_sl, m_sl, p_sl, A, n_sl, B, m_sl, C, n_sl, nr, tol, iwork, dwork, ldwork, info)
+
+
+   !! Create output MATLAB structure
+
+   Ar_mx = mxCreateDoubleMatrix(int(nr, c_size_t), int(nr, c_size_t), mxREAL)
+   Ar(1:nr, 1:nr) => mxGetPr(Ar_mx)
+   Ar = Atmp(1:nr, 1:nr)
+
+   Br_mx = mxCreateDoubleMatrix(int(nr, c_size_t), m, mxREAL)
+   Br(1:nr, 1:m) => mxGetPr(Br_mx)
+   Br = Btmp(1:nr, 1:m)
+
+   Cr_mx = mxCreateDoubleMatrix(p, int(nr, c_size_t), mxREAL)
+   Cr(1:p, 1:nr) => mxGetPr(Cr_mx)
+   Cr = Ctmp(1:p, 1:nr)
+
+   Dr_mx = mxDuplicateArray(D_mx)
+
+   plhs(1) = mxCreateStructMatrix(1_mwSize, 1_mwSize)
+   ignore = mxAddField(plhs(1), "A")
+   call mxSetField(plhs(1), 1_mwIndex, "A", Ar_mx)
+   ignore = mxAddField(plhs(1), "B")
+   call mxSetField(plhs(1), 1_mwIndex, "B", Br_mx)
+   ignore = mxAddField(plhs(1), "C")
+   call mxSetField(plhs(1), 1_mwIndex, "C", Cr_mx)
+   ignore = mxAddField(plhs(1), "D")
+   call mxSetField(plhs(1), 1_mwIndex, "D", Dr_mx)
+end subroutine mexFunction
diff --git a/mex/sources/matlab_mex.F08 b/mex/sources/matlab_mex.F08
index 755b6d5403a20412fb854f9ebcb26d65e7385bed..33f31a4840ee5aeb1c0781c0a88e567cea06e80f 100644
--- a/mex/sources/matlab_mex.F08
+++ b/mex/sources/matlab_mex.F08
@@ -24,7 +24,7 @@
 !   • remove the “value” keywords
 !   • convert input character arrays to character(kind=c_char, len=*)
 
-! Copyright © 2019-2020 Dynare Team
+! Copyright © 2019-2021 Dynare Team
 !
 ! This file is part of Dynare.
 !
@@ -225,6 +225,14 @@ module matlab_mat
      end function mxIsClass_internal
 
      ! Structure
+     type(c_ptr) function mxCreateStructMatrix_internal(m, n, nfields, fieldnames) bind(c, name="mxCreateStructMatrix"//API_VER2)
+       use iso_c_binding
+       import :: mwSize
+       integer(mwSize), intent(in), value :: m, n
+       integer(c_int), intent(in), value :: nfields
+       type(c_ptr), intent(in), value :: fieldnames
+     end function mxCreateStructMatrix_internal
+
      logical(c_bool) function mxIsStruct(pm) bind(c, name="mxIsStruct"//API_VER)
        use iso_c_binding
        type(c_ptr), intent(in), value :: pm
@@ -238,11 +246,25 @@ module matlab_mat
        character(c_char), dimension(*), intent(in) :: fieldname
      end function mxGetField_internal
 
+     subroutine mxSetField_internal(pm, index, fieldname, pvalue) bind(c, name="mxSetField"//API_VER2)
+       use iso_c_binding
+       import :: mwIndex
+       type(c_ptr), intent(in), value :: pm, pvalue
+       integer(mwIndex), intent(in), value :: index
+       character(c_char), dimension(*), intent(in) :: fieldname
+     end subroutine mxSetField_internal
+
      integer(c_int) function mxGetNumberOfFields(pm) bind(c, name="mxGetNumberOfFields"//API_VER)
        use iso_c_binding
        type(c_ptr), intent(in), value :: pm
      end function mxGetNumberOfFields
 
+     integer(c_int) function mxAddField_internal(pm, fieldname) bind(c, name="mxAddField"//API_VER)
+       use iso_c_binding
+       type(c_ptr), intent(in), value :: pm
+       character(c_char), dimension(*), intent(in) :: fieldname
+     end function mxAddField_internal
+
      ! Cell
      logical(c_bool) function mxIsCell(pm) bind(c, name="mxIsCell"//API_VER)
        use iso_c_binding
@@ -326,6 +348,14 @@ contains
     mxIsclass = mxIsclass_internal(pm, classname // c_null_char)
   end function mxIsClass
 
+  ! Do not provide an interface for the last two arguments (initial
+  ! fieldnames), since passing a vector of C character strings is not easy.
+  ! One can simply use mxAddField().
+  type(c_ptr) function mxCreateStructMatrix(m, n)
+    integer(mwSize), intent(in) :: m, n
+    mxCreateStructMatrix = mxCreateStructMatrix_internal(m, n, 0_c_int, c_null_ptr)
+  end function mxCreateStructMatrix
+
   type(c_ptr) function mxGetField(pm, index, fieldname)
     type(c_ptr), intent(in) :: pm
     integer(mwIndex), intent(in) :: index
@@ -333,6 +363,19 @@ contains
     mxGetField = mxGetField_internal(pm, index-1, fieldname // c_null_char)
   end function mxGetField
 
+  subroutine mxSetField(pm, index, fieldname, pvalue)
+    type(c_ptr), intent(in) :: pm, pvalue
+    integer(mwIndex), intent(in) :: index
+    character(kind=c_char, len=*), intent(in) :: fieldname
+    call mxSetField_internal(pm, index-1, fieldname // c_null_char, pvalue)
+  end subroutine mxSetField
+
+  integer(c_int) function mxAddField(pm, fieldname)
+    type(c_ptr), intent(in) :: pm
+    character(kind=c_char, len=*), intent(in) :: fieldname
+    mxAddField = mxAddField_internal(pm, fieldname // c_null_char)
+  end function mxAddField
+
   type(c_ptr) function mxGetCell(pm, index)
     type(c_ptr), intent(in) :: pm
     integer(mwIndex), intent(in) :: index
diff --git a/mex/sources/slicot.F08 b/mex/sources/slicot.F08
new file mode 100644
index 0000000000000000000000000000000000000000..67bc579caa69bb3bb3ccb6c17234b3d8a8e414bf
--- /dev/null
+++ b/mex/sources/slicot.F08
@@ -0,0 +1,42 @@
+! Copyright © 2021 Dynare Team
+!
+! This file is part of Dynare.
+!
+! Dynare is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! Dynare is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+module slicot
+  use iso_fortran_env
+
+#if defined(MATLAB_MEX_FILE) && __SIZEOF_POINTER__ == 8
+  integer, parameter :: slint = int64
+  integer, parameter :: sllog = 8 ! Logical kind, gfortran-specific
+#else
+  integer, parameter :: slint = int32
+  integer, parameter :: sllog = 4 ! Logical kind, gfortran-specific
+#endif
+
+  interface
+     subroutine tb01pd(job, equil, n, m, p, a, lda, b, ldb, c, ldc, &
+          nr, tol, iwork, dwork, ldwork, info)
+       import :: slint, real64
+       character, intent(in) :: job, equil
+       integer(slint), intent(in) :: n, m, p, lda, ldb, ldc, ldwork
+       real(real64), dimension(*), intent(inout) :: a, b, c
+       integer(slint), intent(out) :: nr, info
+       real(real64), intent(in) :: tol
+       integer(slint), dimension(*), intent(out) :: iwork
+       real(real64), dimension(*), intent(out) :: dwork
+     end subroutine tb01pd
+  end interface
+end module slicot