diff --git a/matlab/dr2.m b/matlab/dr2.m
index 662b2c8a9..675513cdc 100644
--- a/matlab/dr2.m
+++ b/matlab/dr2.m
@@ -7,5 +7,5 @@ function ghs2=dr2(ys,dr)
   dr.ys = ys;
   fh = str2func([M_.fname '_static']);
   dr.fbias = 2*feval(fh,dr.ys);
-  dr=dr1(2,dr,0);
+  dr=dr1(dr,0);
   ghs2 = dr.ghs2;
diff --git a/matlab/dynare_estimation.m b/matlab/dynare_estimation.m
index 2bd87a064..8757971a2 100644
--- a/matlab/dynare_estimation.m
+++ b/matlab/dynare_estimation.m
@@ -649,7 +649,7 @@ if (any(bayestopt_.pshape  >0 ) & options_.mh_replic) | (any(bayestopt_.pshape >
   if options_.mh_replic
     metropolis(xparam1,invhess,gend,data,rawdata,bounds);
   end
-  if ~options_.nodiagnostic & options_.mh_replic > 1000 & options_.nblocks > 1
+  if ~options_.nodiagnostic & options_.mh_replic > 1000 & options_.mh_nblck > 1
     McMCDiagnostics;
   end
   %% Here i discard first half of the draws:
diff --git a/matlab/dynare_m.exe b/matlab/dynare_m.exe
index bdf297e2c..bd8fff062 100755
Binary files a/matlab/dynare_m.exe and b/matlab/dynare_m.exe differ
diff --git a/matlab/dynare_resolve.m b/matlab/dynare_resolve.m
index 3791b4e40..03549cf5a 100644
--- a/matlab/dynare_resolve.m
+++ b/matlab/dynare_resolve.m
@@ -1,7 +1,7 @@
 function [A,B,ys,info] = dynare_resolve()
-  global oo_ dr_
+  global oo_
   
-  [dr_,info] = resol(oo_.steady_state,0);
+  [oo_.dr,info] = resol(oo_.steady_state,0);
   
   if info(1) > 0
     A = [];
@@ -10,5 +10,5 @@ function [A,B,ys,info] = dynare_resolve()
     return
   end
   
-  [A,B] = kalman_transition_matrix(dr_);
-  ys = dr_.ys;
\ No newline at end of file
+  [A,B] = kalman_transition_matrix(oo_.dr);
+  ys = oo_.dr.ys;
\ No newline at end of file
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index bbe8dbc5c..1d6f21026 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -1,6 +1,6 @@
 % Copyright (C) 2001 Michel Juillard
 %
-function [x,cheik] = dynare_solve(func,x,varargin)
+function [x,cheik] = dynare_solve(func,x,jacobian_flag,varargin)
   global options_
   
   options_ = set_default_option(options_,'solve_algo',2);
@@ -11,12 +11,16 @@ function [x,cheik] = dynare_solve(func,x,varargin)
       options.MaxFunEvals = 20000;
       options.TolFun=1e-8;
       options.Display = 'off';
-      options.Jacobian = 'on';
+      if jacobian_flag
+	options.Jacobian = 'on';
+      else
+	options.Jacobian = 'on';
+      end
       [x,fval,exitval,output] = fsolve(func,x,options,varargin{:});
       if exitval > 0
-	  	cheik = 0;
+	cheik = 0;
       else
-	  	cheik = 1;
+	cheik = 1;
       end
       return
     else 
@@ -26,14 +30,17 @@ function [x,cheik] = dynare_solve(func,x,varargin)
 
   if options_.solve_algo == 1
     nn = size(x,1) ;
-    [x,cheik]=solve1(func,x,1:nn,1:nn,varargin{:});
+    [x,cheik]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
   elseif options_.solve_algo == 2
     nn = size(x,1) ;
     tolf = eps^(2/3) ;
 
-    %fjac = zeros(nn,nn) ;
-
-    [fvec,fjac] = feval(func,x,varargin{:});
+    if jacobian_flag
+      [fvec,fjac] = feval(func,x,varargin{:});
+    else
+      fvec = feval(func,x,varargin{:});
+      fjac = zeros(nn,nn) ;
+    end
 
     i = find(~isfinite(fvec));
     
@@ -50,17 +57,20 @@ function [x,cheik] = dynare_solve(func,x,varargin)
       return ;
     end
 
-    %dh = max(abs(x),options_.gstep*ones(nn,1))*eps^(1/3);
-    %for j = 1:nn
-    %  xdh = x ;
-    %  xdh(j) = xdh(j)+dh(j) ;
-    %  fjac(:,j) = (feval(func,xdh,varargin{:}) - fvec)./dh(j) ;
-    %end
+    if ~jacobian_flag
+      fjac = zeros(nn,nn) ;
+      dh = max(abs(x),options_.gstep*ones(nn,1))*eps^(1/3);
+      for j = 1:nn
+	xdh = x ;
+	xdh(j) = xdh(j)+dh(j) ;
+	fjac(:,j) = (feval(func,xdh,varargin{:}) - fvec)./dh(j) ;
+      end
+    end
 
     [j1,j2,r,s] = dmperm(fjac);
     
     for i=length(r)-1:-1:1
-      [x,cheik]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),varargin{:});
+      [x,cheik]=solve1(func,x,j1(r(i):r(i+1)-1),j2(r(i):r(i+1)-1),jacobian_flag,varargin{:});
       if cheik
 	error(sprintf('Solve block = %d check = %d\n',i,cheik));
       end
diff --git a/matlab/resol.m b/matlab/resol.m
index 878d38b59..937a78645 100644
--- a/matlab/resol.m
+++ b/matlab/resol.m
@@ -57,7 +57,7 @@ if info(1)
 end
 
 if options_.dr_algo == 1 & options_.order > 1
-  dr.ys = dynare_solve('dr2',ys,dr);
+  dr.ys = dynare_solve('dr2',ys,0,dr);
   dr.fbias = 2*feval([M_.fname '_static'],dr.ys,oo_.exo_steady_state);
   [dr, info1] = dr1(dr,check_flag);
   if info1(1)
diff --git a/matlab/solve1.m b/matlab/solve1.m
index d4f4d749b..29f869612 100644
--- a/matlab/solve1.m
+++ b/matlab/solve1.m
@@ -1,6 +1,6 @@
 % Copyright (C) 2001 Michel Juillard
 %
-function [x,check] = solve1(func,x,j1,j2,varargin)
+function [x,check] = solve1(func,x,j1,j2,jacobian_flag,varargin)
 
   global M_ options_ fjac  
 
@@ -37,17 +37,19 @@ function [x,check] = solve1(func,x,j1,j2,varargin)
   stpmax = stpmx*max([sqrt(x'*x);nn]) ;
   first_time = 1;
   for its = 1:maxit
-
-    %dh = max(abs(x(j2)),options_.gstep*ones(nn,1))*eps^(1/3);
-    
-    %for j = 1:nn
-    %  xdh = x ;
-    %  xdh(j2(j)) = xdh(j2(j))+dh(j) ;
-    %  t = feval(func,xdh,varargin{:});
-    %  fjac(:,j) = (t(j1) - fvec)./dh(j) ;
-    %  g(j) = fvec'*fjac(:,j) ;
-    %end
-    [fvec,fjac] = feval(func,x,varargin{:});
+    if jacobian_flag
+      [fvec,fjac] = feval(func,x,varargin{:});
+    else
+      dh = max(abs(x(j2)),options_.gstep*ones(nn,1))*eps^(1/3);
+      
+      for j = 1:nn
+	xdh = x ;
+	xdh(j2(j)) = xdh(j2(j))+dh(j) ;
+	t = feval(func,xdh,varargin{:});
+	fjac(:,j) = (t(j1) - fvec)./dh(j) ;
+	g(j) = fvec'*fjac(:,j) ;
+      end
+    end
     fvec = fvec(j1);
     fjac = fjac(j1,j2);
     g = (fvec'*fjac)';
diff --git a/matlab/steady_.m b/matlab/steady_.m
index 3ed0d92b9..0ad3e296f 100644
--- a/matlab/steady_.m
+++ b/matlab/steady_.m
@@ -18,7 +18,7 @@ function steady_()
   if exist([M_.fname '_steadystate'])
     [oo_.steady_state,check] = feval([M_.fname '_steadystate'],oo_.steady_state,x);
   else
-    [oo_.staedy_state,check] = dynare_solve([M_.fname '_static'],oo_.steady_state,x);
+    [oo_.steady_state,check] = dynare_solve([M_.fname '_static'],oo_.steady_state,1,x);
   end
 
   if check ~= 0
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 266e67052..83ffd65e0 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -1,7 +1,7 @@
 % Copyright (C) 2001 Michel Juillard
 %
 function info=stoch_simul(var_list)
-global M_ options_ oo_ dr_
+global M_ options_ oo_
 
 global it_
   options_ = set_default_option(options_,'TeX',0);  
@@ -42,7 +42,7 @@ global it_
 
   check_model;
 
-  [dr_, info] = resol(oo_.steady_state,0);
+  [oo_.dr, info] = resol(oo_.steady_state,0);
 
   if info(1)
     print_info(info);
@@ -56,21 +56,21 @@ global it_
     disp(['  Number of variables:         ' int2str(M_.endo_nbr)])
     disp(['  Number of stochastic shocks: ' int2str(M_.exo_nbr)])
     disp(['  Number of state variables:   ' ...
-	  int2str(length(find(dr_.kstate(:,2) <= M_.maximum_lag+1)))])
+	  int2str(length(find(oo_.dr.kstate(:,2) <= M_.maximum_lag+1)))])
     disp(['  Number of jumpers:           ' ...
-	  int2str(length(find(dr_.kstate(:,2) == M_.maximum_lag+2)))])
-    disp(['  Number of static variables:  ' int2str(dr_.nstatic)])
+	  int2str(length(find(oo_.dr.kstate(:,2) == M_.maximum_lag+2)))])
+    disp(['  Number of static variables:  ' int2str(oo_.dr.nstatic)])
     my_title='MATRIX OF COVARIANCE OF EXOGENOUS SHOCKS';
     labels = deblank(M_.exo_names);
     headers = strvcat('Variables',labels);
     lh = size(labels,2)+2;
     table(my_title,headers,labels,M_.Sigma_e,lh,10,6);
     disp(' ')
-    disp_dr(dr_,options_.order,var_list);
+    disp_dr(oo_.dr,options_.order,var_list);
   end
 
   if options_.simul == 0 & options_.nomoments == 0
-    disp_th_moments(dr_,var_list); 
+    disp_th_moments(oo_.dr,var_list); 
   elseif options_.simul == 1
     if options_.periods == 0
       error('STOCH_SIMUL error: number of periods for the simulation isn''t specified')
@@ -80,7 +80,7 @@ global it_
 	    ' than the number of observations to be DROPed'])
       return
     end
-    oo_.endo_simul = simult(repmat(dr_.ys,1,M_.maximum_lag),dr_);
+    oo_.endo_simul = simult(repmat(oo_.dr.ys,1,M_.maximum_lag),oo_.dr);
     dyn2vec;
     if options_.nomoments == 0
       disp_moments(oo_.endo_simul,var_list);
@@ -130,7 +130,7 @@ global it_
     end
     for i=1:M_.exo_nbr
       if SS(i,i) > 1e-13
-	y=irf(dr_,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ...
+	y=irf(oo_.dr,cs(M_.exo_names_orig_ord,i), options_.irf, options_.drop, ...
 	      options_.replic, options_.order);
 	if options_.relative_irf
 	  y = 100*y/cs(i,i); 
@@ -273,28 +273,3 @@ global it_
       fclose(fidTeX);
     end
   end
-  % 02/20/01 MJ oo_.steady_state removed from calling sequence for simult (all in dr_)
-  % 02/23/01 MJ added dyn2vec()
-  % 06/24/01 MJ steady -> steadoo_.endo_simul
-  % 08/28/02 MJ added var_list
-  % 10/09/02 MJ no simulation and theoretical moments for order 1 
-  % 10/14/02 MJ added plot of IRFs
-  % 10/30/02 MJ options_ are now a structure
-  % 01/01/03 MJ added dr_algo
-  % 01/09/03 MJ set default values for options_ (correct absence of autocorr
-  %             when order == 1)
-  % 01/12/03 MJ removed call to steadoo_.endo_simul as already checked in resol()
-  % 02/09/03 MJ oo_.steady_state reset with value declared in initval after computations
-  % 02/18/03 MJ removed above change. oo_.steady_state shouldn't be affected by
-  %             computations in this function
-  %             new option SIMUL computes a stochastic simulation and save
-  %             results in oo_.endo_simul and via dyn2vec
-  % 04/03/03 MJ corrected bug for simulation with M_.maximum_lag > 1
-  % 05/20/03 MJ eliminates exogenous shocks with 0 variance
-  % 05/20/03 MJ don't plot IRF if variation < 1e-10
-  % 11/14/03 MJ corrected bug on number of replications for IRF when
-  %             order=2
-  % 11/22/03 MJ replaced IRFs by orthogonalized IRFs
-  % 08/30/04 SA The maximum number of plots is not constrained for the IRFs and 
-  %			  all the plots are saved in *.eps, *.pdf and *.fig files (added
-  % 09/03/04 SA Tex output for IRFs added
\ No newline at end of file