diff --git a/matlab/dyn_risky_steadystate_solver.m b/matlab/dyn_risky_steadystate_solver.m
index 93c75da4dc4d2b98b92f3f9fafca2668fe0b888e..7fa8a59ffecc9c07644f9f83a805ecbe69ac3617 100644
--- a/matlab/dyn_risky_steadystate_solver.m
+++ b/matlab/dyn_risky_steadystate_solver.m
@@ -62,7 +62,7 @@ function [dr,info] = dyn_risky_steadystate_solver(ys0,M, ...
 %! @end deftypefn
 %@eod:
 
-% Copyright (C) 2001-2018 Dynare Team
+% Copyright (C) 2001-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -165,12 +165,6 @@ if ~isreal(d1) || ~isreal(d2)
     pause
 end
 
-if options.use_dll
-    % In USE_DLL mode, the hessian is in the 3-column sparse representation
-    d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
-                size(d1, 1), size(d1, 2)*size(d1, 2));
-end
-
 if isfield(options,'portfolio') && options.portfolio == 1
     pm = portfolio_model_structure(M,options);
     x = ys(pm.v_p);
@@ -224,12 +218,6 @@ if ~isreal(d1) || ~isreal(d2)
     pause
 end
 
-if options.use_dll
-    % In USE_DLL mode, the hessian is in the 3-column sparse representation
-    d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
-                size(d1, 1), size(d1, 2)*size(d1, 2));
-end
-
 
 gu1 = dr_np.ghu(pm.i_fwrd_g,:);
 
@@ -257,12 +245,6 @@ if ~isreal(d1) || ~isreal(d2)
     pause
 end
 
-if options.use_dll
-    % In USE_DLL mode, the hessian is in the 3-column sparse representation
-    d2 = sparse(d2(:,1), d2(:,2), d2(:,3), ...
-                size(d1, 1), size(d1, 2)*size(d1, 2));
-end
-
 d1_np = d1(pm.eq_np,pm.i_d1_np);
 d2_np = d2(pm.eq_np,pm.i_d2_np);
 
diff --git a/matlab/model_diagnostics.m b/matlab/model_diagnostics.m
index 636713047b9f45f30c708056ef1350b0ee5715c6..2bb2b1990449a96db91df276592b983054c19ed5 100644
--- a/matlab/model_diagnostics.m
+++ b/matlab/model_diagnostics.m
@@ -221,11 +221,6 @@ if ~options.block
                                           exo_simul, ...
                                           M.params, dr.ys, it_);
         end
-        if options.use_dll
-            % In USE_DLL mode, the hessian is in the 3-column sparse representation
-            hessian1 = sparse(hessian1(:,1), hessian1(:,2), hessian1(:,3), ...
-                              size(jacobia_, 1), size(jacobia_, 2)*size(jacobia_, 2));
-        end
     end
 
     if any(any(isinf(jacobia_) | isnan(jacobia_)))
diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m
index 87d88bacf128e8f6cbebe51320344cced1a7d693..85398fa94ce51372471a1cb11204ef90c68d8e11 100644
--- a/matlab/stochastic_solvers.m
+++ b/matlab/stochastic_solvers.m
@@ -128,11 +128,6 @@ elseif local_order == 2
                                       exo_simul, ...
                                       M_.params, dr.ys, it_);
     end
-    if options_.use_dll
-        % In USE_DLL mode, the hessian is in the 3-column sparse representation
-        hessian1 = sparse(hessian1(:,1), hessian1(:,2), hessian1(:,3), ...
-                          size(jacobia_, 1), size(jacobia_, 2)*size(jacobia_, 2));
-    end
     [infrow,infcol]=find(isinf(hessian1));
     if options_.debug
         if ~isempty(infrow)
diff --git a/mex/sources/k_order_perturbation/dynamic_dll.cc b/mex/sources/k_order_perturbation/dynamic_dll.cc
index 4688d66b5bc47d55754b66685195cbc036567e4a..c9a44c6d29d7dd47e798462ecdd8abe85bb29aca 100644
--- a/mex/sources/k_order_perturbation/dynamic_dll.cc
+++ b/mex/sources/k_order_perturbation/dynamic_dll.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2008-2019 Dynare Team
+ * Copyright © 2008-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -19,8 +19,6 @@
 
 #include "dynamic_dll.hh"
 
-#include "dynare_exception.hh"
-
 #include <iostream>
 #include <cassert>
 
@@ -45,34 +43,14 @@ DynamicModelDLL::DynamicModelDLL(const std::string &modName, int ntt_arg, int or
 #endif
                           );
 
-  for (int i = 0; i <= order; i++)
-    {
-      std::string funcname = "dynamic_" + (i == 0 ? "resid" : "g" + std::to_string(i));
-      dynamic_deriv_fct deriv;
-      dynamic_tt_fct tt;
-#if defined(__CYGWIN32__) || defined(_WIN32)
-      deriv = reinterpret_cast<dynamic_deriv_fct>(GetProcAddress(dynamicHinstance, funcname.c_str()));
-      tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, (funcname + "_tt").c_str()));
-#else
-      deriv = reinterpret_cast<dynamic_deriv_fct>(dlsym(dynamicHinstance, funcname.c_str()));
-      tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, (funcname + "_tt").c_str()));
-#endif
-      if (!deriv || !tt)
-        {
-#if defined(__CYGWIN32__) || defined(_WIN32)
-          FreeLibrary(dynamicHinstance);
-#else
-          dlclose(dynamicHinstance);
-#endif
-          throw DynareException(__FILE__, __LINE__, "Error when loading symbols from " + fName
-#if !defined(__CYGWIN32__) && !defined(_WIN32)
-                                + ": " + dlerror()
-#endif
-                                );
-        }
-      dynamic_deriv.push_back(deriv);
-      dynamic_tt.push_back(tt);
-    }
+  dynamic_tt.resize(order+1);
+
+  std::tie(dynamic_resid, dynamic_tt[0]) = getSymbolsFromDLL<dynamic_resid_or_g1_fct>("dynamic_resid", fName);
+  std::tie(dynamic_g1, dynamic_tt[1]) = getSymbolsFromDLL<dynamic_resid_or_g1_fct>("dynamic_g1", fName);
+
+  dynamic_higher_deriv.resize(std::max(0, order-1));
+  for (int i = 2; i <= order; i++)
+    std::tie(dynamic_higher_deriv[i-2], dynamic_tt[i]) = getSymbolsFromDLL<dynamic_higher_deriv_fct>("dynamic_g" + std::to_string(i), fName);
 
   tt = std::make_unique<double[]>(ntt);
 }
@@ -95,11 +73,16 @@ void
 DynamicModelDLL::eval(const Vector &y, const Vector &x, const Vector &modParams, const Vector &ySteady,
                       Vector &residual, std::vector<TwoDMatrix> &md) noexcept(false)
 {
-  assert(md.size() == dynamic_deriv.size()-1);
+  assert(md.size() == dynamic_tt.size()-1);
 
-  for (size_t i = 0; i < dynamic_deriv.size(); i++)
+  for (size_t i = 0; i < dynamic_tt.size(); i++)
     {
       dynamic_tt[i](y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get());
-      dynamic_deriv[i](y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), i == 0 ? residual.base() : md[i-1].base());
+      if (i == 0)
+        dynamic_resid(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), residual.base());
+      else if (i == 1)
+        dynamic_g1(y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), md[0].base());
+      else
+        dynamic_higher_deriv[i-2](y.base(), x.base(), 1, modParams.base(), ySteady.base(), 0, tt.get(), &md[i-1].get(0, 0), &md[i-1].get(0, 1), &md[i-1].get(0, 2));
     }
 }
diff --git a/mex/sources/k_order_perturbation/dynamic_dll.hh b/mex/sources/k_order_perturbation/dynamic_dll.hh
index 693e8b2b9d9a7aa583496fe5d25baed02bdafbf0..e9d9aa8a6ff9c3afb7242139944db24ee8432b15 100644
--- a/mex/sources/k_order_perturbation/dynamic_dll.hh
+++ b/mex/sources/k_order_perturbation/dynamic_dll.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright © 2008-2019 Dynare Team
+ * Copyright © 2008-2020 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -31,11 +31,15 @@
 
 #include <string>
 #include <memory>
+#include <utility>
+
+#include "dynare_exception.hh"
 
 #include "dynamic_abstract_class.hh"
 
 using dynamic_tt_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, double *T);
-using dynamic_deriv_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *deriv);
+using dynamic_resid_or_g1_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *resid_or_g1);
+using dynamic_higher_deriv_fct = void (*)(const double *y, const double *x, int nb_row_x, const double *params, const double *steady_state, int it_, const double *T, double *g_i, double *g_j, double *g_v);
 
 /**
  * creates pointer to Dynamic function inside <model>_dynamic.dll
@@ -45,7 +49,8 @@ class DynamicModelDLL : public DynamicModelAC
 {
 private:
   std::vector<dynamic_tt_fct> dynamic_tt;
-  std::vector<dynamic_deriv_fct> dynamic_deriv;
+  dynamic_resid_or_g1_fct dynamic_resid, dynamic_g1;
+  std::vector<dynamic_higher_deriv_fct> dynamic_higher_deriv; // Index 0 is g2
 #if defined(_WIN32) || defined(__CYGWIN32__)
   HINSTANCE dynamicHinstance; // DLL instance pointer in Windows
 #else
@@ -53,6 +58,34 @@ private:
 #endif
   std::unique_ptr<double[]> tt; // Vector of temporary terms
 
+  template<typename T>
+  std::pair<T, dynamic_tt_fct>
+  getSymbolsFromDLL(const std::string &funcname, const std::string &fName)
+  {
+    dynamic_tt_fct tt;
+    T deriv;
+#if defined(__CYGWIN32__) || defined(_WIN32)
+    deriv = reinterpret_cast<T>(GetProcAddress(dynamicHinstance, funcname.c_str()));
+    tt = reinterpret_cast<dynamic_tt_fct>(GetProcAddress(dynamicHinstance, (funcname + "_tt").c_str()));
+#else
+      deriv = reinterpret_cast<T>(dlsym(dynamicHinstance, funcname.c_str()));
+      tt = reinterpret_cast<dynamic_tt_fct>(dlsym(dynamicHinstance, (funcname + "_tt").c_str()));
+#endif
+      if (!deriv || !tt)
+        {
+#if defined(__CYGWIN32__) || defined(_WIN32)
+          FreeLibrary(dynamicHinstance);
+#else
+          dlclose(dynamicHinstance);
+#endif
+          throw DynareException(__FILE__, __LINE__, "Error when loading symbols from " + fName
+#if !defined(__CYGWIN32__) && !defined(_WIN32)
+                                + ": " + dlerror()
+#endif
+                                );
+        }
+      return { deriv, tt };
+  }
 public:
   // construct and load Dynamic model DLL
   explicit DynamicModelDLL(const std::string &fname, int ntt_arg, int order);
diff --git a/preprocessor b/preprocessor
index ad5e196d3064fd3b431b0a0dd34c66f42b52dbc6..c4351166a9828d67e555c617fb97dfc00a1bc880 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit ad5e196d3064fd3b431b0a0dd34c66f42b52dbc6
+Subproject commit c4351166a9828d67e555c617fb97dfc00a1bc880