diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 11310e984ea3d0a8dc9b6011b429915d67c7814f..921a1d53e5bc8f85de69b7d88d992d5ce14d2462 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -105,7 +105,7 @@ elseif options_.solve_algo == 2 || options_.solve_algo == 4
 
     if ~jacobian_flag
         fjac = zeros(nn,nn) ;
-        dh = max(abs(x),options_.gstep*ones(nn,1))*eps^(1/3);
+        dh = max(abs(x),options_.gstep(1)*ones(nn,1))*eps^(1/3);
         for j = 1:nn
             xdh = x ;
             xdh(j) = xdh(j)+dh(j) ;
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 406cee2514d2d94677cbf7bd2e0b04f095ad6bed..0d50155c41a2840a61e2c9f62cc01daded3e87ed 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -45,7 +45,9 @@ options_.dynatol.x = 1e-5;
 options_.maxit_ = 10;
 options_.slowc = 1;
 options_.timing = 0;
-options_.gstep = 1e-2;
+options_.gstep = ones(2,1);
+options_.gstep(1) = 1e-2;
+options_.gstep(2) = 1.0;
 options_.scalv = 1;
 options_.debug = 0;
 options_.initval_file = 0;
diff --git a/matlab/hessian.m b/matlab/hessian.m
index 3b8fdb8b9c119b7f6e3906090dbd1f86df311d51..9e1825cd54b92857250867eaadd3076a186fe16b 100644
--- a/matlab/hessian.m
+++ b/matlab/hessian.m
@@ -39,7 +39,7 @@ if ~isa(func, 'function_handle')
     func = str2func(func);
 end
 n=size(x,1);
-h1=max(abs(x),sqrt(gstep)*ones(n,1))*eps^(1/6);
+h1=max(abs(x),sqrt(gstep(1))*ones(n,1))*eps^(1/6)*gstep(2);
 h_1=h1;
 xh1=x+h1;
 h1=xh1-x;
diff --git a/matlab/solve1.m b/matlab/solve1.m
index 6b3dcec813e1618faadb71d0569de76888e4c0bf..104ccf3328ed86ff9ec4e486299ed3dac8a5a96f 100644
--- a/matlab/solve1.m
+++ b/matlab/solve1.m
@@ -77,7 +77,7 @@ for its = 1:maxit
         fvec = fvec(j1);
         fjac = fjac(j1,j2);
     else
-        dh = max(abs(x(j2)),options_.gstep*ones(nn,1))*eps^(1/3);
+        dh = max(abs(x(j2)),options_.gstep(1)*ones(nn,1))*eps^(1/3);
         
         for j = 1:nn
             xdh = x ;