diff --git a/matlab/csolve.m b/matlab/csolve.m
index 5f95fe128..e6400b830 100644
--- a/matlab/csolve.m
+++ b/matlab/csolve.m
@@ -43,19 +43,23 @@ af0=sum(abs(f0));
 af00=af0;
 itct=0;
 while ~done
-   disp([af00-af0 crit*max(1,af0)])
+%   disp([af00-af0 crit*max(1,af0)])
    if itct>3 & af00-af0<crit*max(1,af0) & rem(itct,2)==1
       randomize=1;
    else
       if ~analyticg
-         if isempty(varargin)
-            grad = (feval(FUN,x*ones(1,nv)+tvec)-f0*ones(1,nv))/delta;
-         else
-            grad = (feval(FUN,x*ones(1,nv)+tvec,varargin{:})-f0*ones(1,nv))/delta;
-         end
+% $$$          if isempty(varargin)
+% $$$             grad = (feval(FUN,x*ones(1,nv)+tvec)-f0*ones(1,nv))/delta;
+% $$$          else
+% $$$             grad = (feval(FUN,x*ones(1,nv)+tvec,varargin{:})-f0*ones(1,nv))/delta;
+% $$$          end
+	grad = zeros(nv,nv);
+	for i=1:nv
+	  grad(:,i) = (feval(FUN,x+tvec(:,i),varargin{:})-f0)/delta;
+	end
       else % use analytic gradient
 	   %         grad=feval(gradfun,x,varargin{:});
-	   [f0,grad] = feval(FUN,x,varargin{:});
+	   [f0,grad] = feval(gradfun,x,varargin{:});
       end
       if isreal(grad)
          if rcond(grad)<1e-12
diff --git a/matlab/dynare_solve.m b/matlab/dynare_solve.m
index 5a18e8923..b63dd8c81 100644
--- a/matlab/dynare_solve.m
+++ b/matlab/dynare_solve.m
@@ -8,9 +8,10 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
   if options_.solve_algo == 0
     if ~isempty(which('fsolve')) & sscanf(version('-release'),'%d') >= 13;
       options=optimset('fsolve');
-      options.MaxFunEvals = 20000;
+      options.MaxFunEvals = 50000;
+      options.MaxIter = 2000;
       options.TolFun=1e-8;
-      options.Display = 'off';
+      options.Display = 'iter';
       if jacobian_flag
 	options.Jacobian = 'on';
       else
@@ -51,7 +52,7 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
       error('exiting ...')
     end
     
-    f = 0.5*fvec'*fvec ;
+%    f = 0.5*fvec'*fvec ;
 
     if max(abs(fvec)) < tolf
       return ;
@@ -80,7 +81,11 @@ function [x,info] = dynare_solve(func,x,jacobian_flag,varargin)
       [x,info]=solve1(func,x,1:nn,1:nn,jacobian_flag,varargin{:});
     end
   elseif options_.solve_algo == 3
-      [x,info] = csolve(func,x,'grad_ss',1e-6,500,varargin{:});
+    if jacobian_flag
+      [x,info] = csolve(func,x,func,1e-6,500,varargin{:});
+    else
+      [x,info] = csolve(func,x,[],1e-6,500,varargin{:});
+    end 
   end
 %    fvec1 = feval(func,x,varargin{:})