diff --git a/matlab/numgrad2_.m b/matlab/numgrad2_.m
new file mode 100644
index 000000000..3140a2ff7
--- /dev/null
+++ b/matlab/numgrad2_.m
@@ -0,0 +1,60 @@
+function [g, badg] = numgrad2_(fcn,f0,x,epsilon,scale,varargin)
+% function [g badg] = numgrad2(fcn,xvarargin)
+
+% Original file downloaded from:
+% http://sims.princeton.edu/yftp/optimize/mfiles/numgrad.m
+
+% Copyright (C) 1993-2007 Christopher Sims
+% Copyright (C) 2012 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see .
+
+fh = NaN;
+
+delta = epsilon*max(abs(x));
+n = length(x);
+g = zeros(n,1);
+
+badg=0;
+goog=1;
+g0 = 0;
+
+for i=1:n
+ xiold = x(i);
+ h = step_length_correction(xiold,scale,i)*delta;
+ x(i) = xiold + h;
+ [fh,junk1,junk2,cost_flag] = feval(fcn, x, varargin{:});
+ if cost_flag
+ g0 = (fh - f0)/h;
+ else
+ x(i) = xiold - h;
+ [fh,junk1,junk2,cost_flag] = feval(fcn, x, varargin{:});
+ if cost_flag
+ g0 = (f0-fh)/h;
+ else
+ goog = 0;
+ end
+ end
+ if goog && abs(g0)< 1e15
+ g(i) = g0;
+ else
+ disp('bad gradient ------------------------')
+ % fprintf('Gradient w.r.t. %3d: %10g\n',i,g0)
+ g(i) = 0;
+ badg = 1;
+ end
+ x(i) = xiold;
+end
\ No newline at end of file
diff --git a/matlab/numgrad3_.m b/matlab/numgrad3_.m
new file mode 100644
index 000000000..67f126036
--- /dev/null
+++ b/matlab/numgrad3_.m
@@ -0,0 +1,72 @@
+function [g, badg, f0, f1, f2] = numgrad3_(fcn,f0,x,epsilon,scale,varargin)
+% Computes the gradient of the objective function fcn using a three points
+% formula if possible.
+%
+% Adapted from Sims' numgrad routine.
+%
+% See section 25.3.4 in Abramovitz and Stegun (1972, Tenth Printing, December) Handbook of Mathematical Functions.
+% http://www.math.sfu.ca/~cbm/aands/
+
+% Original file downloaded from:
+% http://sims.princeton.edu/yftp/optimize/mfiles/numgrad.m
+
+% Copyright (C) 1993-2007 Christopher Sims
+% Copyright (C) 2008-2012 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see .
+
+f1 = NaN;
+f2 = NaN;
+
+delta = epsilon;
+n = length(x);
+g = zeros(n,1);
+
+badg=0;
+goog=1;
+
+zgrad = 1;
+
+for i=1:n
+ xiold = x(i);
+ h = step_length_correction(xiold,scale,i)*delta;
+ x(i) = xiold + h;
+ [f1,junk1,junk2,cost_flag1] = feval(fcn, x, varargin{:});
+ x(i) = xiold - h;
+ [f2,junk1,junk2,cost_flag2] = feval(fcn, x, varargin{:});
+ if cost_flag1 && cost_flag2
+ g0 = (f1 - f2) / (2*h);
+ if zgrad && f1>f0 && f2>f0 % Note that this condition is consistent with a minimization problem!
+ g0 = 0;
+ end
+ else
+ if cost_flag1
+ g0 = (f1-f0)/h;
+ elseif cost_flag2
+ g0 = (f0-f2)/h;
+ else
+ goog=0;
+ end
+ end
+ if goog && abs(g0)< 1e15
+ g(i)=g0;
+ else
+ disp('bad gradient ------------------------')
+ g(i)=0;
+ badg=1;
+ end
+ x(i) = xiold;
+end
\ No newline at end of file
diff --git a/matlab/numgrad5_.m b/matlab/numgrad5_.m
new file mode 100644
index 000000000..ac9e45655
--- /dev/null
+++ b/matlab/numgrad5_.m
@@ -0,0 +1,91 @@
+function [g, badg, f0, f1, f2, f3, f4] = numgrad5_(fcn,f0,x,epsilon,scale,varargin)
+% Computes the gradient of the objective function fcn using a five points
+% formula if possible.
+%
+% Adapted from Sims' numgrad.m routine.
+%
+% See section 25.3.6 Abramovitz and Stegun (1972, Tenth Printing, December) Handbook of Mathematical Functions.
+% http://www.math.sfu.ca/~cbm/aands/
+%
+% TODO Try Four points formula when cost_flag3=0 or cost_flag4=0.
+
+% Original file downloaded from:
+% http://sims.princeton.edu/yftp/optimize/mfiles/numgrad.m
+
+% Copyright (C) 1993-2007 Christopher Sims
+% Copyright (C) 2008-2012 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare. If not, see .
+
+f1 = NaN;
+f2 = NaN;
+f3 = NaN;
+f4 = NaN;
+
+delta = epsilon;
+n=length(x);
+tvec=delta*eye(n);
+g=zeros(n,1);
+
+badg=0;
+goog=1;
+
+zgrad = 1;
+
+for i=1:n
+ xiold = x(i);
+ h = step_length_correction(xiold,scale,i)*delta;
+ x(i) = xiold+h;
+ [f1,junk1,junk2,cost_flag1] = feval(fcn, x, varargin{:});
+ x(i) = xiold-h;
+ [f2,junk1,junk2,cost_flag2] = feval(fcn, x, varargin{:});
+ if ~cost_flag1 || ~cost_flag2
+ cost_flag3 = 0;
+ cost_flag4 = 0;
+ disp('numgrad:: I cannot use the five points formula!!')
+ else
+ x(i) = xiold+2*h;
+ [f3,junk1,junk2,cost_flag3] = feval(fcn, x, varargin{:});
+ x(i) = xiold-2*h;
+ [f4,junk1,junk2,cost_flag4] = feval(fcn, x, varargin{:});
+ end
+ if cost_flag1 && cost_flag2 && cost_flag3 && cost_flag4% Five Points formula
+ g0 = (8*(f1 - f2)+ f4-f3) / (12*h);
+ if zgrad && f0