diff --git a/matlab/trust_region.m b/matlab/trust_region.m
deleted file mode 100644
index 1edb87732c000b535ed78715082555aba0d73640..0000000000000000000000000000000000000000
--- a/matlab/trust_region.m
+++ /dev/null
@@ -1,245 +0,0 @@
-function [x,check] = trust_region(fcn,x0,j1,j2,jacobian_flag,gstep,tolf,tolx,maxiter,debug,varargin)
-% Solves systems of non linear equations of several variables, using a
-% trust-region method.
-%
-% INPUTS
-% fcn: name of the function to be solved
-% x0: guess values
-% j1: equations index for which the model is solved
-% j2: unknown variables index
-% jacobian_flag=1: jacobian given by the 'func' function
-% jacobian_flag=0: jacobian obtained numerically
-% gstep increment multiplier in numercial derivative
-% computation
-% tolf tolerance for residuals
-% tolx tolerance for solution variation
-% maxiter maximum number of iterations
-% debug debug flag
-% varargin: list of arguments following bad_cond_flag
-%
-% OUTPUTS
-% x: results
-% check=1: the model can not be solved
-%
-% SPECIAL REQUIREMENTS
-% none
-
-% Copyright (C) 2008-2012 VZLU Prague, a.s.
-% Copyright (C) 2014 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 <http://www.gnu.org/licenses/>.
-%
-% Initial author: Jaroslav Hajek <highegg@gmail.com>, for GNU Octave
-
-if (ischar (fcn))
- fcn = str2func (fcn);
-end
-
-n = length(j1);
-
-% These defaults are rather stringent. I think that normally, user
-% prefers accuracy to performance.
-
-macheps = eps (class (x0));
-
-niter = 1;
-
-x = x0;
-info = 0;
-
-% Initial evaluation.
-% Handle arbitrary shapes of x and f and remember them.
-fvec = fcn (x, varargin{:});
-fvec = fvec(j1);
-fn = norm (fvec);
-jcn = nan(n, 1);
-
-% Outer loop.
-while (niter < maxiter && ~info)
-
- % Calculate function value and Jacobian (possibly via FD).
- if jacobian_flag
- [fvec, fjac] = fcn (x, varargin{:});
- fvec = fvec(j1);
- fjac = fjac(j1,j2);
- else
- dh = max(abs(x(j2)),gstep(1)*ones(n,1))*eps^(1/3);
-
- for j = 1:n
- xdh = x ;
- xdh(j2(j)) = xdh(j2(j))+dh(j) ;
- t = fcn(xdh,varargin{:});
- fjac(:,j) = (t(j1) - fvec)./dh(j) ;
- g(j) = fvec'*fjac(:,j) ;
- end
- end
-
- % Get column norms, use them as scaling factors.
- for j = 1:n
- jcn(j) = norm(fjac(:,j));
- end
- if (niter == 1)
- dg = jcn;
- dg(dg == 0) = 1;
- else
- % Rescale adaptively.
- % FIXME: the original minpack used the following rescaling strategy:
- % dg = max (dg, jcn);
- % but it seems not good if we start with a bad guess yielding Jacobian
- % columns with large norms that later decrease, because the corresponding
- % variable will still be overscaled. So instead, we only give the old
- % scaling a small momentum, but do not honor it.
-
- dg = max (0.1*dg, jcn);
- end
-
- if (niter == 1)
- xn = norm (dg .* x(j2));
- % FIXME: something better?
- delta = max (xn, 1);
- end
-
- % Get trust-region model (dogleg) minimizer.
- s = - dogleg (fjac, fvec, dg, delta);
- w = fvec + fjac * s;
-
- sn = norm (dg .* s);
- if (niter == 1)
- delta = min (delta, sn);
- end
-
- x2 = x;
- x2(j2) = x2(j2) + s;
- fvec1 = fcn (x2, varargin{:});
- fvec1 = fvec1(j1);
- fn1 = norm (fvec1);
-
- if (fn1 < fn)
- % Scaled actual reduction.
- actred = 1 - (fn1/fn)^2;
- else
- actred = -1;
- end
-
- % Scaled predicted reduction, and ratio.
- t = norm (w);
- if (t < fn)
- prered = 1 - (t/fn)^2;
- ratio = actred / prered;
- else
- prered = 0;
- ratio = 0;
- end
-
- % Update delta.
- if (ratio < 0.1)
- delta = 0.5*delta;
- if (delta <= 1e1*macheps*xn)
- % Trust region became uselessly small.
- info = -3;
- break;
- end
- elseif (abs (1-ratio) <= 0.1)
- delta = 1.4142*sn;
- elseif (ratio >= 0.5)
- delta = max (delta, 1.4142*sn);
- end
-
- if (ratio >= 1e-4)
- % Successful iteration.
- x(j2) = x(j2) + s;
- xn = norm (dg .* x(j2));
- fvec = fvec1;
- fn = fn1;
- end
-
- niter = niter + 1;
-
-
- % Tests for termination conditions. A mysterious place, anything
- % can happen if you change something here...
-
- % The rule of thumb (which I'm not sure M*b is quite following)
- % is that for a tolerance that depends on scaling, only 0 makes
- % sense as a default value. But 0 usually means uselessly long
- % iterations, so we need scaling-independent tolerances wherever
- % possible.
-
- % FIXME -- why tolf*n*xn? If abs (e) ~ abs(x) * eps is a vector
- % of perturbations of x, then norm (fjac*e) <= eps*n*xn, i.e. by
- % tolf ~ eps we demand as much accuracy as we can expect.
- disp([niter fn ratio])
- if (fn <= tolf*n*xn)
- info = 1;
- % The following tests done only after successful step.
- elseif (ratio >= 1e-4)
- % This one is classic. Note that we use scaled variables again,
- % but compare to scaled step, so nothing bad.
- if (sn <= tolx*xn)
- info = 2;
- % Again a classic one. It seems weird to use the same tolf
- % for two different tests, but that's what M*b manual appears
- % to say.
- elseif (actred < tolf)
- info = 3;
- end
- end
-end
-
-check = ~info;
-end
-
-
-% Solve the double dogleg trust-region least-squares problem:
-% Minimize norm(r*x-b) subject to the constraint norm(d.*x) <= delta,
-% x being a convex combination of the gauss-newton and scaled gradient.
-
-% TODO: error checks
-% TODO: handle singularity, or leave it up to mldivide?
-
-function x = dogleg (r, b, d, delta)
-% Get Gauss-Newton direction.
-x = r \ b;
-xn = norm (d .* x);
-if (xn > delta)
- % GN is too big, get scaled gradient.
- s = (r' * b) ./ d;
- sn = norm (s);
- if (sn > 0)
- % Normalize and rescale.
- s = (s / sn) ./ d;
- % Get the line minimizer in s direction.
- tn = norm (r*s);
- snm = (sn / tn) / tn;
- if (snm < delta)
- % Get the dogleg path minimizer.
- bn = norm (b);
- dxn = delta/xn; snmd = snm/delta;
- t = (bn/sn) * (bn/xn) * snmd;
- t = t - dxn * snmd^2 + sqrt ((t-dxn)^2 + (1-dxn^2)*(1-snmd^2));
- alpha = dxn*(1-snmd^2) / t;
- else
- alpha = 0;
- end
- else
- alpha = delta / xn;
- snm = 0;
- end
- % Form the appropriate convex combination.
- x = alpha * x + ((1-alpha) * min (snm, delta)) * s;
-end
-end
-