diff --git a/matlab/PosteriorIRF_core1.m b/matlab/PosteriorIRF_core1.m
index f264cf1a7..002d77ba8 100644
--- a/matlab/PosteriorIRF_core1.m
+++ b/matlab/PosteriorIRF_core1.m
@@ -194,7 +194,7 @@ while fpar<B
     end
     if MAX_nirfs_dsgevar
         IRUN = IRUN+1;
-        [fval,info,cost_flag,junk1,junk2,SteadyState,junk3,PHI,SIGMAu,iXX] =  dsge_var_likelihood_1(deep',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
+        [fval,junk1,junk2,cost_flag,SteadyState,junk3,info,PHI,SIGMAu,iXX] =  dsge_var_likelihood(deep',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,bounds,oo_);
         dsge_prior_weight = M_.params(strmatch('dsge_prior_weight',M_.param_names));
         DSGE_PRIOR_WEIGHT = floor(dataset_.nobs*(1+dsge_prior_weight));
         SIGMA_inv_upper_chol = chol(inv(SIGMAu*dataset_.nobs*(dsge_prior_weight+1)));
diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m
index dbd1c6db4..642023e0a 100644
--- a/matlab/initial_estimation_checks.m
+++ b/matlab/initial_estimation_checks.m
@@ -117,9 +117,9 @@ ana_deriv = DynareOptions.analytic_derivation;
 DynareOptions.analytic_derivation=0;
 if ~isequal(DynareOptions.mode_compute,11) || ...
         (isequal(DynareOptions.mode_compute,11) && isequal(DynareOptions.order,1))
-  [fval,a,b,c,d] = objective_function(xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults);
+  [fval,junk1,junk2,a,b,c,d] = feval(objective_function,xparam1,DynareDataset,DatasetInfo,DynareOptions,Model,EstimatedParameters,BayesInfo,BoundsInfo,DynareResults);
 else 
-    a=0;
+    b=0;
     fval = 0;
 end
 if DynareOptions.debug
@@ -127,10 +127,10 @@ if DynareOptions.debug
 end
 DynareOptions.analytic_derivation=ana_deriv;
 
-if strcmp(func2str(objective_function),'non_linear_dsge_likelihood')
-    info = d;
+if DynareOptions.dsge_var || strcmp(func2str(objective_function),'non_linear_dsge_likelihood')
+    info = b;
 else
-    info = a;
+    info = d;
 end
 
 % if DynareOptions.mode_compute==5
diff --git a/matlab/optimization/numgrad5_.m b/matlab/optimization/numgrad5_.m
index 1d03e62ac..f696f4c8c 100644
--- a/matlab/optimization/numgrad5_.m
+++ b/matlab/optimization/numgrad5_.m
@@ -1,4 +1,4 @@
-function [g, badg] = numgrad5_(fcn,f0,x,penalty,epsilon,varargin)
+function [g, badg] = numgrad5(fcn,f0,x,epsilon,varargin)
 % Computes the gradient of the objective function fcn using a five points
 % formula if possible.
 %
@@ -48,13 +48,13 @@ for i=1:n
     xiold = x(i);
     h = step_length_correction(xiold,scale,i)*delta;
     x(i) = xiold+h;
-    [f1,junk1,junk2,cost_flag1] = penalty_objective_function(x, fcn, penalty, varargin{:});
+    [f1,junk1,junk2,cost_flag1] = feval(fcn, x, penalty, varargin{:});
     x(i) = xiold-h;
-    [f2,junk1,junk2,cost_flag2] = penalty_objective_function(x, fcn, penalty, varargin{:});
+    [f2,junk1,junk2,cost_flag2] = feval(fcn, x, penalty, varargin{:});
     x(i) = xiold+2*h;
-    [f3,junk1,junk2,cost_flag3] = penalty_objective_function(x, fcn, penalty, varargin{:});
+    [f3,junk1,junk2,cost_flag3] = feval(fcn, x, penalty, varargin{:});
     x(i) = xiold-2*h;
-    [f4,junk1,junk2,cost_flag4] = penalty_objective_function(x, fcn, penalty, varargin{:});
+    [f4,junk1,junk2,cost_flag4] = feval(fcn, x, penalty, varargin{:});
     if f0<f1 && f1<f3 && f0<f2 && f2<f4
         g0 = 0;
     else