diff --git a/matlab/mr_gstep.m b/matlab/mr_gstep.m
index 06e6965ae..ce6448fbe 100644
--- a/matlab/mr_gstep.m
+++ b/matlab/mr_gstep.m
@@ -1,5 +1,5 @@
-function [f0, x, ig] = mr_gstep(func0,x,htol0,varargin)
-% function [f0, x] = mr_gstep(func0,x,htol0,varargin)
+function [f0, x, ig] = mr_gstep(init,x,func0,htol0,varargin)
+% function [f0, x] = mr_gstep(init,x,func0,htol0,varargin)
 % 
 % Gibbs type step in optimisation
 
@@ -23,21 +23,22 @@ function [f0, x, ig] = mr_gstep(func0,x,htol0,varargin)
 global bayestopt_ options_
 persistent h1 
 
-gstep_ = options_.gstep;
-if nargin<3, 
+n=size(x,1);
+
+if init,
+    gstep_ = options_.gstep;
+    h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/4);
+    return
+end
+if nargin<4, 
     htol = 1.e-6;
 else
     htol = htol0;
 end
 func = str2func(func0);
 f0=feval(func,x,varargin{:});
-n=size(x,1);
 h2=bayestopt_.ub-bayestopt_.lb;
 
-if isempty(h1),
-    h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/4);
-end
-
 xh1=x;
 f1=zeros(size(f0,1),n);
 f_1=f1;
diff --git a/matlab/mr_hessian.m b/matlab/mr_hessian.m
index 18e9f6fd8..efa5490c0 100644
--- a/matlab/mr_hessian.m
+++ b/matlab/mr_hessian.m
@@ -1,5 +1,5 @@
-function [hessian_mat, gg, htol1, ihh, hh_mat0] = mr_hessian(func,x,hflag,htol0,varargin)
-%  [hessian_mat, gg, htol1, ihh, hh_mat0] = mr_hessian(func,x,hflag,htol0,varargin)
+function [hessian_mat, gg, htol1, ihh, hh_mat0] = mr_hessian(init,x,func,hflag,htol0,varargin)
+%  [hessian_mat, gg, htol1, ihh, hh_mat0] = mr_hessian(init,x,func,hflag,htol0,varargin)
 %
 %  numerical gradient and Hessian, with 'automatic' check of numerical
 %  error 
@@ -44,19 +44,20 @@ function [hessian_mat, gg, htol1, ihh, hh_mat0] = mr_hessian(func,x,hflag,htol0,
 global options_ bayestopt_
 persistent h1 htol
 
-gstep_=options_.gstep;
-if isempty(htol), htol = 1.e-4; end
+n=size(x,1);
+if init,
+    gstep_=options_.gstep;
+    htol = 1.e-4; 
+%h1=max(abs(x),gstep_*ones(n,1))*eps^(1/3);
+%h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/6);
+    h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/4);
+    return,
+end
 func = str2func(func);
 [f0, ff0]=feval(func,x,varargin{:});
-n=size(x,1);
 h2=bayestopt_.ub-bayestopt_.lb;
 hmax=bayestopt_.ub-x;
 hmax=min(hmax,x-bayestopt_.lb);
-%h1=max(abs(x),gstep_*ones(n,1))*eps^(1/3);
-%h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/6);
-if isempty(h1),
-    h1=max(abs(x),sqrt(gstep_)*ones(n,1))*eps^(1/4);
-end
 
 h1 = min(h1,0.5.*hmax);
 
diff --git a/matlab/newrat.m b/matlab/newrat.m
index 3a256e279..0e11f12e4 100644
--- a/matlab/newrat.m
+++ b/matlab/newrat.m
@@ -59,8 +59,12 @@ func_hh = [func0,'_hh'];
 func = str2func(func0);
 fval0=feval(func,x,varargin{:});
 fval=fval0;
+% initialize mr_gstep and mr_hessian
+mr_gstep(1,x);
+mr_hessian(1,x);
+
 if isempty(hh)
-    [dum, gg, htol0, igg, hhg]=mr_hessian(func_hh,x,flagit,htol,varargin{:});
+    [dum, gg, htol0, igg, hhg]=mr_hessian(0,x,func_hh,flagit,htol,varargin{:});
     hh0 = reshape(dum,nx,nx);
     hh=hhg;
     if min(eig(hh0))<0,
@@ -121,7 +125,7 @@ while norm(gg)>gtol & check==0 & jit<nit,
             iggx(find(ig),find(ig)) = inv( hhx(find(ig),find(ig)) );
             [fvala x0 fc retcode] = csminit(func0,x0,fval,ggx,0,iggx,varargin{:});
         end
-        [fvala, x0, ig] = mr_gstep(func0,x0,htol,varargin{:});
+        [fvala, x0, ig] = mr_gstep(0,x0,func0,htol,varargin{:});
         nig=[nig ig];
         if (fval-fvala)<gibbstol*(fval0(icount)-fval),
             igibbs=0;
@@ -162,7 +166,7 @@ while norm(gg)>gtol & check==0 & jit<nit,
         if flagit==2,
             hh=hh0;
         elseif flagg>0,
-            [dum, gg, htol0, igg, hhg]=mr_hessian(func_hh,xparam1,flagg,ftol0,varargin{:});   
+            [dum, gg, htol0, igg, hhg]=mr_hessian(0,xparam1,func_hh,flagg,ftol0,varargin{:});   
             if flagg==2,
                 hh = reshape(dum,nx,nx);
                 ee=eig(hh);
@@ -202,7 +206,7 @@ while norm(gg)>gtol & check==0 & jit<nit,
             catch
                 save m1 x fval0 nig 
             end
-            [dum, gg, htol0, igg, hhg]=mr_hessian(func_hh,xparam1,flagit,htol,varargin{:});
+            [dum, gg, htol0, igg, hhg]=mr_hessian(0,xparam1,func_hh,flagit,htol,varargin{:});
             if htol0>htol, %ftol,
                            %ftol=htol0;
                 htol=htol0;