diff --git a/matlab/method_of_moments/method_of_moments.m b/matlab/method_of_moments/method_of_moments.m
index c83bae21bedf054b18cc2f0b05649ce6befca048..337b8e0bcc733339d9a4cc5b61dcce35b21a1afa 100644
--- a/matlab/method_of_moments/method_of_moments.m
+++ b/matlab/method_of_moments/method_of_moments.m
@@ -93,7 +93,6 @@ function [oo_, options_mom_, M_] = method_of_moments(bayestopt_, options_, oo_,
 % - [ ] deal with measurement errors (once @wmutschl has implemented this in identification toolbox)
 % - [ ] improve check for duplicate moments by using the cellfun and unique functions
 % - [ ] dirname option to save output to different directory not yet implemented
-% - [ ] add analytic_jacobian option for mode_compute 4 and 101
 
 % The TeX option crashes MATLAB R2014a run with "-nodisplay" option
 % (as is done from the testsuite).
@@ -352,7 +351,7 @@ options_mom_.vector_output= false;           % specifies whether the objective f
 
 optimizer_vec=[options_mom_.mode_compute;num2cell(options_mom_.additional_optimizer_steps)]; % at each stage one can possibly use different optimizers sequentially
 
-analytic_jacobian_optimizers = [1, 3, 13]; %these are currently supported, see to-do list
+analytic_jacobian_optimizers = [1, 3, 4, 13, 101]; %these are currently supported, see to-do list
 
 % -------------------------------------------------------------------------
 % Step 1d: Other options that need to be initialized
@@ -914,24 +913,11 @@ for stage_iter=1:size(options_mom_.mom.weighting_matrix,1)
             end
             if strcmp(options_mom_.mom.mom_method,'GMM') && options_mom_.mom.analytic_jacobian && ismember(optimizer_vec{optim_iter},analytic_jacobian_optimizers) %do this only for gradient-based optimizers
                 options_mom_.mom.compute_derivs = true;
-                objective_function1 = str2func('method_of_moments_objective_function_gradient_helper');
-                switch optimizer_vec{optim_iter}
-                    case 1
-                        options_mom_.optim_opt = [optim_opt0, ',''GradObj'',''on'''];   % make sure GradObj option is on for fmincon
-                    case 3
-                        options_mom_.optim_opt = [optim_opt0, ',''GradObj'',''on'''];   % make sure GradObj option is on for fmincon
-                    case 13
-                        options_mom_.optim_opt = [optim_opt0, ',''Jacobian'',''on'''];  % make sure Jacobian option is on for lsqnonlin
-                end
-                if strcmp(options_mom_.optim_opt(1),',')
-                    options_mom_.optim_opt(1) = []; %remove the comma if optim_opt was empty
-                end
             else
                 options_mom_.mom.compute_derivs = false;
-                objective_function1 = objective_function;
             end
             
-            [xparam1, fval, exitflag] = dynare_minimize_objective(objective_function1, xparam0, optimizer_vec{optim_iter}, options_mom_, [Bounds.lb Bounds.ub], bayestopt_laplace.name, bayestopt_laplace, [],...
+            [xparam1, fval, exitflag] = dynare_minimize_objective(objective_function, xparam0, optimizer_vec{optim_iter}, options_mom_, [Bounds.lb Bounds.ub], bayestopt_laplace.name, bayestopt_laplace, [],...
                                                                   Bounds, oo_, estim_params_, M_, options_mom_);
             if options_mom_.vector_output
                 fval = fval'*fval;
diff --git a/matlab/method_of_moments/method_of_moments_objective_function.m b/matlab/method_of_moments/method_of_moments_objective_function.m
index 3a8548e76105de126ddeb8fdb2980fc20cccd2bf..9aecc87a25363097946a8f89681d8442c674fa5f 100644
--- a/matlab/method_of_moments/method_of_moments_objective_function.m
+++ b/matlab/method_of_moments/method_of_moments_objective_function.m
@@ -1,5 +1,5 @@
-function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_, df] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
-% [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
+function [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
+% [fval, info, exit_flag, df, junk1, oo_, M_, options_mom_] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
 % -------------------------------------------------------------------------
 % This function evaluates the objective function for GMM/SMM estimation
 % =========================================================================
@@ -15,14 +15,13 @@ function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_, df] = meth
 %   o fval:                     value of the quadratic form of the moment difference (except for lsqnonlin, where this is done implicitly)
 %   o info:                     vector storing error code and penalty 
 %   o exit_flag:                0 if error, 1 if no error
-%   o junk1:                    empty matrix required for optimizer interface
-%   o junk2:                    empty matrix required for optimizer interface
+%   o df:                       analytical parameter Jacobian of the quadratic form of the moment difference (for GMM only)
+%   o junk1:                    empty matrix required for optimizer interface (Hessian would go here)
 %   o oo_:                      structure containing the results with the following updated fields:
 %      - mom.model_moments       [numMom x 1] vector with model moments
 %      - mom.Q                   value of the quadratic form of the moment difference
 %   o M_:                       Matlab's structure describing the model
 %   o options_mom_:             structure information about all settings (specified by the user, preprocessor, and taken from global options_)
-%   o df:                       analytical parameter Jacobian of the quadratic form of the moment difference (for GMM only)
 % -------------------------------------------------------------------------
 % This function is called by
 %  o method_of_moments.m
@@ -59,14 +58,18 @@ function [fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_, df] = meth
 %------------------------------------------------------------------------------
 % 0. Initialization of the returned variables and others...
 %------------------------------------------------------------------------------
-if options_mom_.vector_output == 1
-    if options_mom_.mom.penalized_estimator
-        df           = nan(size(oo_.mom.data_moments,1)+length(xparam1),length(xparam1));
+if options_mom_.mom.compute_derivs && options_mom_.mom.analytic_jacobian
+    if options_mom_.vector_output == 1
+        if options_mom_.mom.penalized_estimator
+            df           = nan(size(oo_.mom.data_moments,1)+length(xparam1),length(xparam1));
+        else
+            df           = nan(size(oo_.mom.data_moments,1),length(xparam1));
+        end
     else
-        df           = nan(size(oo_.mom.data_moments,1),length(xparam1));
+        df           = nan(1,length(xparam1));
     end
 else
-    df           = nan(1,length(xparam1));
+    df=[]; %required to be empty by e.g. newrat
 end
 junk1        = [];
 junk2        = [];
diff --git a/matlab/method_of_moments/method_of_moments_objective_function_gradient_helper.m b/matlab/method_of_moments/method_of_moments_objective_function_gradient_helper.m
deleted file mode 100644
index f1c88af499ae2d8367c2a62ddd9a65b3a96230cb..0000000000000000000000000000000000000000
--- a/matlab/method_of_moments/method_of_moments_objective_function_gradient_helper.m
+++ /dev/null
@@ -1,55 +0,0 @@
-function [out1, out2] = method_of_moments_objective_function_gradient_helper(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
-% [out1, out2] = method_of_moments_objective_function_gradient_helper(xparam1, Bounds, oo_, estim_params_, M_, options_mom_)
-% -------------------------------------------------------------------------
-% This helper function evaluates the objective function for GMM/SMM estimation and
-% outputs the function value fval at first and the gradient df at second place,
-% needed for gradient-based optimizers if analytic_jacobian option is set
-% =========================================================================
-% INPUTS
-%   o xparam1:                  current value of estimated parameters as returned by set_prior()
-%   o Bounds:                   structure containing parameter bounds
-%   o oo_:                      structure for results
-%   o estim_params_:            structure describing the estimated_parameters
-%   o M_                        structure describing the model
-%   o options_mom_:             structure information about all settings (specified by the user, preprocessor, and taken from global options_)
-% -------------------------------------------------------------------------
-% OUTPUTS: dependent on the optimizer calling this function, see below
-% -------------------------------------------------------------------------
-% This function calls
-%  o method_of_moments_objective_function.m
-% =========================================================================
-% Copyright (C) 2020-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/>.
-% -------------------------------------------------------------------------
-% This function is called by
-%  o method_of_moments.m
-%  o dynare_minimize_objective.m
-% -------------------------------------------------------------------------
-% This function calls
-%  o method_of_moments_objective_function
-% -------------------------------------------------------------------------
-% Author(s): 
-% o Willi Mutschler (willi@mutschler.eu)
-% =========================================================================
-[fval, info, exit_flag, junk1, junk2, oo_, M_, options_mom_, df] = method_of_moments_objective_function(xparam1, Bounds, oo_, estim_params_, M_, options_mom_);
-
-% mode_compute=1|3|13 require the gradient as second output argument
-out1=fval;
-out2=df;
-
-end%main function end
-
diff --git a/matlab/optimization/analytic_gradient_wrapper.m b/matlab/optimization/analytic_gradient_wrapper.m
new file mode 100644
index 0000000000000000000000000000000000000000..a74484f9964c0e819daa7b3f696d039a6ba232c1
--- /dev/null
+++ b/matlab/optimization/analytic_gradient_wrapper.m
@@ -0,0 +1,34 @@
+function [fval, grad, hess, exit_flag]=analytic_gradient_wrapper(x, fcn, varargin)
+%function [fval, grad, hess, exitflag]=analytic_gradient_wrapper(x, fcn, varargin)
+% Encapsulates an objective function to be minimized for use with Matlab
+% optimizers
+%
+% INPUTS
+% - x             [double]    n*1 vector of instrument values.
+% - fcn           [fhandle]   objective function.
+% - varagin       [cell]      additional parameters for fcn.
+%
+% OUTPUTS
+% - fval          [double]    scalar, value of the objective function at x.
+% - grad                      gradient of the objective function
+% - hess                      Hessian of the objective function
+% - exit_flag     [integer]   scalar, flag returned by
+
+% Copyright (C) 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/>.
+
+[fval, info, exit_flag, grad, hess] = fcn(x, varargin{:});
\ No newline at end of file