diff --git a/matlab/csminwel.m b/matlab/csminwel.m
index 846863edd..03ff36384 100644
--- a/matlab/csminwel.m
+++ b/matlab/csminwel.m
@@ -1,4 +1,4 @@
-function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit,method,varargin)
+function [fh,xh,gh,H,itct,fcount,retcodeh] = csminwel(fcn,x0,H0,grad,crit,nit,method,epsilon,varargin)
 %[fhat,xhat,ghat,Hhat,itct,fcount,retcodehat] = csminwel(fcn,x0,H0,grad,crit,nit,method,varargin)
 % fcn:   string naming the objective function to be minimized
 % x0:    initial value of the parameter vector
@@ -72,13 +72,13 @@ end
 if NumGrad
     switch method
       case 2
-        [g,badg] = numgrad(fcn, f0, x0, varargin{:});
+        [g,badg] = numgrad(fcn, f0, x0, epsilon, varargin{:});
         badg
         g
       case 3
-        [g,badg] = numgrad3(fcn, f0, x0, varargin{:});
+        [g,badg] = numgrad3(fcn, f0, x0, epsilon, varargin{:});
       case 5
-        [g,badg] = numgrad5(fcn, f0, x0, varargin{:});    
+        [g,badg] = numgrad5(fcn, f0, x0, epsilon, varargin{:});    
     end
 else
     [g,badg] = feval(grad,x0,varargin{:});
@@ -121,11 +121,11 @@ while ~done
             if NumGrad
                 switch method 
                   case 2
-                    [g1 badg1] = numgrad(fcn, f1, x1, varargin{:});
+                    [g1 badg1] = numgrad(fcn, f1, x1, epsilon, varargin{:});
                   case 3
-                    [g1 badg1] = numgrad3(fcn, f1, x1, varargin{:});
+                    [g1 badg1] = numgrad3(fcn, f1, x1, epsilon, varargin{:});
                   case 5
-                    [g1,badg1] = numgrad5(fcn, f1, x1, varargin{:});             
+                    [g1,badg1] = numgrad5(fcn, f1, x1, epsilon, varargin{:});             
                 end
             else
                 [g1 badg1] = feval(grad,x1,varargin{:});
@@ -155,11 +155,11 @@ while ~done
                     if NumGrad
                         switch method
                           case 2
-                            [g2 badg2] = numgrad(fcn, f2, x2, varargin{:});
+                            [g2 badg2] = numgrad(fcn, f2, x2, epsilon, varargin{:});
                           case 3
-                            [g2 badg2] = numgrad3(fcn, f2, x2, varargin{:});
+                            [g2 badg2] = numgrad3(fcn, f2, x2, epsilon, varargin{:});
                           case 5
-                            [g2,badg2] = numgrad5(fcn, f2, x2, varargin{:});                   
+                            [g2,badg2] = numgrad5(fcn, f2, x2, epsilon, varargin{:});                   
                         end
                     else
                         [g2 badg2] = feval(grad,x2,varargin{:});
@@ -190,11 +190,11 @@ while ~done
                             if NumGrad
                                 switch method
                                   case 2
-                                    [g3 badg3] = numgrad(fcn, f3, x3, varargin{:});
+                                    [g3 badg3] = numgrad(fcn, f3, x3, epsilon, varargin{:});
                                   case 3
-                                    [g3 badg3] = numgrad3(fcn, f3, x3, varargin{:});
+                                    [g3 badg3] = numgrad3(fcn, f3, x3, epsilon, varargin{:});
                                   case 5
-                                    [g3,badg3] = numgrad5(fcn, f3, x3, varargin{:});                         
+                                    [g3,badg3] = numgrad5(fcn, f3, x3, epsilon, varargin{:});                         
                                 end
                             else
                                 [g3 badg3] = feval(grad,x3,varargin{:});
@@ -254,11 +254,11 @@ while ~done
             if NumGrad
                 switch method
                   case 2
-                    [gh,badgh] = numgrad(fcn, fh, xh, varargin{:});
+                    [gh,badgh] = numgrad(fcn, fh, xh, epsilon, varargin{:});
                   case 3
-                    [gh,badgh] = numgrad3(fcn, fh, xh, varargin{:});
+                    [gh,badgh] = numgrad3(fcn, fh, xh, epsilon, varargin{:});
                   case 5
-                    [gh,badgh] = numgrad5(fcn, fh, xh, varargin{:});
+                    [gh,badgh] = numgrad5(fcn, fh, xh, epsilon, varargin{:});
                 end
             else
                 [gh badgh] = feval(grad, xh,varargin{:});
diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index de8ac28b7..4a4111f89 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -12,7 +12,7 @@ function dynare_estimation_1(var_list_,dname)
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright (C) 2003-2009 Dynare Team
+% Copyright (C) 2003-2009, 2010 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -400,12 +400,12 @@ if options_.mode_compute > 0 & options_.posterior_mode_estimation
         verbose = 2;
         if ~options_.bvar_dsge
             [fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ...
-                csminwel('DsgeLikelihood',xparam1,H0,[],crit,nit,options_.gradient_method,gend,data,data_index,number_of_observations,no_more_missing_observations);
+                csminwel('DsgeLikelihood',xparam1,H0,[],crit,nit,options_.gradient_method,options_.gradient_epsilon,gend,data,data_index,number_of_observations,no_more_missing_observations);
             disp(sprintf('Objective function at mode: %f',fval))
             disp(sprintf('Objective function at mode: %f',DsgeLikelihood(xparam1,gend,data,data_index,number_of_observations,no_more_missing_observations)))
         else
             [fval,xparam1,grad,hessian_csminwel,itct,fcount,retcodehat] = ...
-                csminwel('DsgeVarLikelihood',xparam1,H0,[],crit,nit,options_.gradient_method,gend);
+                csminwel('DsgeVarLikelihood',xparam1,H0,[],crit,nit,options_.gradient_method,options_.gradient_epsilon,gend);
             disp(sprintf('Objective function at mode: %f',fval))
             disp(sprintf('Objective function at mode: %f',DsgeVarLikelihood(xparam1,gend)))
         end
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 70914a0fe..fc2282edd 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -182,6 +182,7 @@ options_.subdraws = [];
 options_.unit_root_vars = [];
 options_.use_mh_covariance_matrix = 0;
 options_.gradient_method = 2;
+options_.gradient_epsilon = 1e-6;
 options_.posterior_sampling_method = 'random_walk_metropolis_hastings';
 options_.proposal_distribution = 'rand_multivariate_normal';
 options_.student_degrees_of_freedom = 3;
diff --git a/matlab/maximize_prior_density.m b/matlab/maximize_prior_density.m
index ca0587683..3ff454399 100644
--- a/matlab/maximize_prior_density.m
+++ b/matlab/maximize_prior_density.m
@@ -15,7 +15,7 @@ function [xparams,lpd,hessian] = ...
 %   lpd           [double]  scalar, value of the logged prior density at the mode.
 %   hessian       [double]  matrix, Hessian matrix at the prior mode.
 
-% Copyright (C) 2009 Dynare Team
+% Copyright (C) 2009, 2010 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -40,7 +40,7 @@ verbose = 2;
 gradient_method = 2;
 
 [lpd,xparams,grad,hessian,itct,fcount,retcodehat] = ...
-    csminwel('minus_logged_prior_density',iparams,H0,[],crit,nit,gradient_method, ... 
+    csminwel('minus_logged_prior_density',iparams,H0,[],crit,nit,gradient_method, options_.gradient_epsilon, ... 
              prior_shape, prior_hyperparameter_1, prior_hyperparameter_2, prior_inf_bound, prior_sup_bound);
 
 lpd = -lpd;
diff --git a/matlab/numgrad.m b/matlab/numgrad.m
index 13c2f329c..69809efa3 100644
--- a/matlab/numgrad.m
+++ b/matlab/numgrad.m
@@ -1,4 +1,4 @@
-function [g, badg] = numgrad(fcn,f0,x,varargin)
+function [g, badg] = numgrad(fcn,f0,x,epsilon,varargin)
 % function [g badg] = numgrad(fcn,xvarargin)
 
 % Original file downloaded from:
diff --git a/matlab/numgrad3.m b/matlab/numgrad3.m
index 6128e1984..b62d0ac85 100644
--- a/matlab/numgrad3.m
+++ b/matlab/numgrad3.m
@@ -1,4 +1,4 @@
-1function [g, badg, f0, f1, f2] = numgrad3(fcn,f0,x,varargin)
+1function [g, badg, f0, f1, f2] = numgrad3(fcn,f0,x,epsilon,varargin)
 % Computes the gradient of the objective function fcn using a three points
 % formula if possible.
 %
diff --git a/matlab/numgrad5.m b/matlab/numgrad5.m
index 1125826e3..c1bfb07ec 100644
--- a/matlab/numgrad5.m
+++ b/matlab/numgrad5.m
@@ -1,4 +1,4 @@
-function [g, badg, f0, f1, f2, f3, f4] = numgrad5(fcn,f0,x,varargin)
+function [g, badg, f0, f1, f2, f3, f4] = numgrad5(fcn,f0,x,epsilon,varargin)
 % Computes the gradient of the objective function fcn using a five points
 % formula if possible.
 %
diff --git a/matlab/osr1.m b/matlab/osr1.m
index 29fbe2ab6..6574a67e2 100644
--- a/matlab/osr1.m
+++ b/matlab/osr1.m
@@ -1,6 +1,6 @@
 function osr1(i_params,i_var,weights)
 
-% Copyright (C) 2005-2009 Dynare Team
+% Copyright (C) 2005-2009, 2010 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -95,7 +95,7 @@ crit = 1e-7;
 nit = 1000;
 verbose = 2;
 
-[f,p]=csminwel('osr_obj',t0,H0,[],crit,nit,options_.gradient_method,i_params,...
+[f,p]=csminwel('osr_obj',t0,H0,[],crit,nit,options_.gradient_method,options_.gradient_epsilon,i_params,...
                inv_order_var(i_var),weights(i_var,i_var));
 
 %  options = optimset('fminunc');
diff --git a/matlab/simulated_moments_estimation.m b/matlab/simulated_moments_estimation.m
index e1a5b58c1..b4f41eaaa 100644
--- a/matlab/simulated_moments_estimation.m
+++ b/matlab/simulated_moments_estimation.m
@@ -149,11 +149,11 @@ vb = 2;
 % Minimization of the objective function.
 if nargin==3
     [fval,param,grad,hessian_csminwel,itct,fcount,retcodehat] = ...
-        csminwel('smm_objective',xparam,H0,[],ct,it,2,sample_moments,weighting_matrix,options);
+        csminwel('smm_objective',xparam,H0,[],ct,it,2,options_.gradient_epsilon,sample_moments,weighting_matrix,options);
 
 elseif nargin>3
     [fval,param,grad,hessian_csminwel,itct,fcount,retcodehat] = ...
-        csminwel('smm_objective',xparam,H0,[],ct,it,2,sample_moments,weighting_matrix,options,parallel);
+        csminwel('smm_objective',xparam,H0,[],ct,it,2,options_.gradient_epsilon,sample_moments,weighting_matrix,options,parallel);
 end